SUBROUTINE simulatemun()

  USE Commonvars
  USE Random
  USE Statistics
  USE Nelder_Mead

  IMPLICIT NONE
  INCLUDE 'mpif.h'
  
  REAL(8), ALLOCATABLE :: simout(:,:,:), simoutt(:,:,:,:), simouttt(:,:,:,:), betart(:), betanort(:), beta_i(:), beta_c(:), var_base(:,:,:), fr_stolen_terms_t(:,:,:,:), fr_stolen_terms(:,:,:,:)
  !, betapmayor_base(:,:,:,:,:,:,:,:,:,:,:), betat_base(:)
  INTEGER, ALLOCATABLE :: terms_t(:,:,:), terms(:,:,:)
  REAL(8), EXTERNAL    :: zeroin, wealth_equivalent
  INTEGER, PARAMETER   :: N_other_var = 8
  REAL(8)              :: msavingf, ch_saving, stealf, qconsf, mwagef,  mnewsavingf, newstealf, newpconsf, newqconsf, yy, fundsf, utilf, ts1, ts2, ts3, zzpuf, zzprf, uud(6), uud_age(2), pauditt(1:2), &
                          emaxr, emaxnor, numt(Nmomtemp), dent(Nmomtemp), num(Nmomtemp), den(Nmomtemp), other_var_numt(N_other_var), other_var_dent(N_other_var), &
                          other_var_num(N_other_var), other_var_den(N_other_var), other_var(N_other_var), fr_stolen, dum_steal, dum_steal_x, vemaxpf, &
                          ln_relqcons, d_run, fine_fun, rundec_netW, zzpr_shock_t, run_shock_t, thres_frac_stolen, newstealf_nome, &
                          ability, potential_wage, wage_pastmayor, msaving_bf, med_wealth(3), frac_fine_p, frac_fine_c, will_to_pay, av_will_to_pay_t(Nactmun,Nsimterms), &
                          av_will_to_pay(Nactmun,Nsimterms), av_fr_stolen_t(Nactmun,Nsimterms), av_fr_stolen(Nactmun,Nsimterms), av_frst_terms(Nnterms), av_frst, &
                          ELNqcons, temp_var(Nvar), uselect, Pwin(2), emax_voter_i, emax_voter_c, alpha_nstt, alpha_eat, d_past, Pr_win_i, &
                          simmom_p(Nmomtemp,1), popsizet, pp, Phi_arg, arg_exp, pnorm, d_audit, av_qcons, diff_qcons, p50_frstolen(1), p75_qcons(1), ptest(1), &
                          fr_stolent(Nactmun), qconst(Nactmun), av_sav_2t, av_log_ab_2t, std_log_ab_2t, second_mom_log_ab_2t, &
                          prob_age_2t(Retage-1), prob_edu_2t(Neduc), prob_mun_2t(Nmun), multi_slope, V_slope, notmayor_base_value, ELNqcons_base, vemax_ini, &
                          vemax_base, vemax_policy, wealth_star, a, b, fa, fb, fr_wealth, av_fr_stolen_small, av_fr_stolen_large, av_wtp, av_wtp_resources, &
                          av_wtp_income, years_to_death, pc_cost_prau_t(Nmun), pc_cost_wage_t(Nmun), prob_conv_vec(2), meas_err_shock, &
                          mean_lnstealing, mean2_lnstealing, pr_ch_age_vec(2)
  INTEGER              :: i, j, ii, jj, ieduv, index_steal_x, tt, outelect, outfine, auditf,auditf3, thief3, mayorage, mayorage_elec, ch_age, incumb_age, incumb_educ, nterms, ipastf, Totterms, &
                          aenoc_mun, aud, aud3, ipau, ipauf, totloop, ntasks, restasks, counter, iw, is, ifu, ipr, ist, imun, ieducf, ch_educ, isav, infeasible, dft, &
                          nsimt, thief, ub_sav, ub_sav_ch, ub_savi, ub_sav_med, med_edu(3), med_age(3), med_ipast(3), med_edu_i, med_ipast_i, i_num, i_den, voter_age, &
                          elap, elapv, count_av, count_low_edu, count_age(Retage-1), i1, i2, i3, iselect, outmun, nprocs_t, myrank_t, iboot, out_ch_age
  CHARACTER(len=100)   :: sim_var_names(Nmomtemp), data_var_names(Nmomtemp)
  EXTERNAL             :: objfcn_ml, notmayordecision, rmultinomial, mayordecision, pnorm

  INTERFACE
     SUBROUTINE sfun(n, x, f, g)
       IMPLICIT NONE
       INTEGER, PARAMETER      :: dp = SELECTED_REAL_KIND(12, 60)
       INTEGER, INTENT(IN)     :: n
       REAL (dp), INTENT(IN)   :: x(:)
       REAL (dp), INTENT(OUT)  :: f, g(:)
     END SUBROUTINE sfun
     SUBROUTINE sfun1(n, x, f, g)
       IMPLICIT NONE
       INTEGER, PARAMETER      :: dp = SELECTED_REAL_KIND(12, 60)
       INTEGER, INTENT(IN)     :: n
       REAL (dp), INTENT(IN)   :: x(:)
       REAL (dp), INTENT(OUT)  :: f, g(:)
     END SUBROUTINE sfun1
  END INTERFACE


  flag_simulation = 1

  CALL cpu_time(ts1)
  
  IF (flag_policy_evaluation == 0 .OR. (flag_policy_evaluation == 1 .AND. flag_optimal_auditprob == 1)) THEN
     nprocs_t = nprocs
     myrank_t = myrank
  ELSE
     nprocs_t = nprocs_new
     myrank_t = my_new_comm_rank
  END IF

!!$  IF (myrank_t == 0 .AND. flag_policy_evaluation == 1) write(*,'(a,1i4,11f18.14)') 'param inputs', nboot, theta, delta_fine, alpha_q, cost_run_frac, upower, alpha_ie, alpha_ea, eta, alpha_st, pappeal(1), sigma_steal

  !!$  nprocs_t = 1

  !Compute the size of all loop
  totloop = Nsim
  ntasks = INT(totloop/nprocs_t)
  restasks = totloop - ntasks*nprocs_t

  ALLOCATE(simout(Nactmun,Nsimterms,Nvar))
  simout = 0d0

  IF (flag_simulate_data == 1 .OR. (flag_policy_evaluation == 1 .AND. nboot==1)) THEN
     ALLOCATE(simoutt(Nactmun,Nsim,Nsimterms,Nvar))
     simoutt = 0d0
  END IF

  ALLOCATE(terms_t(Nactmun,Nsim,Nsimterms), terms(Nactmun,Nsim,Nsimterms), fr_stolen_terms_t(Nnterms,Nactmun,Nsim,Nsimterms), fr_stolen_terms(Nnterms,Nactmun,Nsim,Nsimterms))
  terms_t = 0
  terms = 0
  fr_stolen_terms_t = 0d0
  fr_stolen_terms = 0d0
  
  av_will_to_pay_t = 0d0
  av_fr_stolen_t = 0d0
  !In 2019 in Brazil, median age is 33.2 (2=32-37 in our code) and average education is 9.3 (1=high school), We use median age, education, and wealth conditional on the size of the municipality, computed in structest_Feb_2020.do
  !We allow the median value to vary with municipality size
  med_edu(1:2) = 1
  med_edu(3) = 1
  med_age(1) = 2
  med_age(2:3) = 2
  med_ipast(:) = 1

  ! cost of the audit policy in wealth terms
  IF (cost_policy == 1) THEN
     IF (flag_prau_policy == 1 .OR. flag_optimal_prau_policy == 1 .OR. flag_term_optimal_prau_policy == 1 .OR. flag_norun_optimal_prau_policy == 1 .OR. flag_optimal_wage_policy == 1 .OR. flag_term_norun_optimal_prau_policy == 1 .OR. flag_wage_optimal_prau_policy == 1) THEN
        DO imun = 1, Nmun
           pc_cost_prau_t(imun) = DBLE(Tend - med_age(imun) + 1)*pc_cost_prau
!!$           IF (myrank == 0) write(*,'(a,3f20.10)') 'pc_cost_prau_t(imun)', pc_cost_prau_t(imun), DBLE(Tend - med_age(imun) + 1), pc_cost_prau
        END DO
     END IF
     IF (flag_wage_policy == 1 .OR. flag_optimal_wage_policy == 1 .OR. flag_wage_optimal_prau_policy == 1) THEN
        DO imun = 1, Nmun
           pc_cost_wage_t(imun) = DBLE(Tend - med_age(imun) + 1)*pc_cost_wage(imun)
!!$           IF (myrank == 0) write(*,'(a,3f20.10)') 'pc_cost_wage_t(imun)', pc_cost_wage_t(imun), DBLE(Tend - med_age(imun) + 1), pc_cost_wage(imun)
        END DO
     END IF
  END IF
  
!!$  Median and mean wealth in Brazil:
!!$  The numbers are in US dollars, we have to multiply by X to get number in Brazilian Reals; What exchange rate do we use? The current one is $1 = 5.2335 Reals, the average for 2019 is $1 = 3.9457 Reals
!!$  for 2001-2010 (mean) for 2010 (median): https://www.credit-suisse.com/media/assets/corporate/docs/about-us/research/publications/credit-suisse-global-wealth-databook.pdf
!!$  For 2019: https://www.credit-suisse.com/media/assets/corporate/docs/about-us/research/publications/global-wealth-databook-2019.pdf
!!$  for 2019: https://en.wikipedia.org/wiki/List_of_countries_by_wealth_per_adult
!!$  med_wealth(1) = 218.248d0
!!$  med_wealth(2) = 287.3192d0
!!$  med_wealth(3) = 430.7853d0
  med_wealth(1) = 47.38653233d0 ! = 14,242 * 0.633892551 * 5.2489 / norm (average per-adult brazilian wealth in US dollars in 2006 in Table 2.4 times price index to deflate to 2000 times the average exchange rate in April 1st 2020) !66.9727943d0 ! = 24,289 * 0.525315827 * 5.2489 / norm (average per-adult brazilian wealth in US dollars in 2010 in Table 2.4 times price index to deflate to 2000 times the average exchange rate in April 1st 2020) ! 92.921235d0 ! = 23,550 * 3.9457/norm (average per-adult brazilian wealth in US dollars in 2019 in Table 21 times the average exchange rate in 2019).
  med_wealth(2) = 47.38653233d0 !66.9727943d0 ! = 24,289 * 0.525315827 * 5.2489 / norm (average per-adult brazilian wealth in US dollars in 2010 in Table 2.4 times price index to deflate to 2000 times the average exchange rate in April 1st 2020).
  med_wealth(3) = 47.38653233d0 !66.9727943d0
!!$  !Given the choosen wealth, education, and age, total initial resources in the computation of wtp = 156.272929 = 92.921235d0*1.05+pmwage*12
!!$  !Wtp as fraction of wealth should therefore be computed as wtp/156.272929
!!$  !We read in the residents' values for the base case
  
  !Find the point on the saving domain where we need to compute the approximated value function
  DO isav = 1, Nwealth
     IF (med_wealth(2) <  wealth(isav)) THEN
        ub_sav_med = isav
        EXIT
     END IF
     ub_sav_med = Nwealth
  END DO
  
  Totterms = Nsimterms
  numt = 0d0
  dent = 0d0
  num = 0d0
  den = 0d0
  
  other_var_numt = 0d0
  other_var_dent = 0d0
  other_var_num = 0d0
  other_var_den = 0d0

  frac_fine_c = 0d0

!!$        Nsimt = 1
  Nsimt = nsim

  ! The following integers are only used in the computation of emax and not in the simulation,
  ! but we need them in mayordecision.f90
  iw = 0
  is = 0
  ifu = 0
  ipr = 0

  ! If we control for observable heterogeneity at the mayor level
  IF (flag_no_selec_sav == 1 .OR. flag_no_selec_age == 1 .OR. flag_no_selec_educ == 1 .OR. flag_no_selec_ab == 1 .OR. flag_no_selec_funds == 1 .OR. flag_no_selec_zzpr == 1 .OR. flag_no_selec_mun == 1) THEN

!!$     ALLOCATE(var_base(Nsim,Nactmun,Nsimterms,2), sav_inc(Nsim*Nactmun,Nmun))
     ALLOCATE(var_base(Nsim*Nactmun,Nmun,6))
     var_base = 0d0
     
     OPEN(unit=18, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_base_Jan2023.dat")  

     av_sav_2t = 0d0
     av_log_ab_2t = 0d0
     second_mom_log_ab_2t = 0d0
     count_av = 0
     count_low_edu = 0
     count_age = 0
     i1 = 0
     i2 = 0
     i3 = 0
     DO j = 1, Nsim
        DO i = 1, Nactmun
           DO tt = 1, Nsimterms
              
              READ(18,*) temp_var
            
              IF (temp_var(15) > 0.5d0 .AND. tt ==2) THEN
                 av_sav_2t = av_sav_2t + temp_var(12)
                 av_log_ab_2t = av_log_ab_2t + LOG(temp_var(14))
                 second_mom_log_ab_2t = second_mom_log_ab_2t + LOG(temp_var(14))**2d0
                 count_av = count_av + 1
                 IF (temp_var(13) < 1.5d0) count_low_edu =  count_low_edu + 1
                 DO ii = 1, Retage - 1
                    IF (ABS(temp_var(27)-DBLE(ii)) < 0.1d0) count_age(ii) = count_age(ii)+1
                 END DO

                 imun = INT(temp_var(23))
                 IF (imun == 1) THEN
                    i1 = i1 + 1
                    ! Savings
                    var_base(i1,imun,1) = temp_var(12)
                    ! Ability
                    var_base(i1,imun,2) = temp_var(14)
                    ! Age
                    var_base(i1,imun,3) = temp_var(27)
                    ! Education
                    var_base(i1,imun,4) = temp_var(13)
                    ! Private Inputs
                    var_base(i1,imun,5) = temp_var(10)
                    ! Funds
                    var_base(i1,imun,6) = temp_var(9)                    
!!$                    sav_inc(i1,imun) =  temp_var(12)
                 ELSE IF (imun == 2) THEN
                    i2 = i2 + 1
                    ! Savings
                    var_base(i2,imun,1) = temp_var(12)
                    ! Ability
                    var_base(i2,imun,2) = temp_var(14)
                    ! Age
                    var_base(i2,imun,3) = temp_var(27)
                    ! Education
                    var_base(i2,imun,4) = temp_var(13)
                    ! Private Inputs
                    var_base(i2,imun,5) = temp_var(10)
                    ! Funds
                    var_base(i2,imun,6) = temp_var(9)
!!$                    sav_inc(i2,imun) =  temp_var(12)
                 ELSE IF (imun == 3) THEN
                    i3 = i3 + 1
                    ! Savings
                    var_base(i3,imun,1) = temp_var(12)
                    ! Ability
                    var_base(i3,imun,2) = temp_var(14)
                    ! Age
                    var_base(i3,imun,3) = temp_var(27)
                    ! Education
                    var_base(i3,imun,4) = temp_var(13)
                    ! Private Inputs
                    var_base(i3,imun,5) = temp_var(10)
                    ! Funds
                    var_base(i3,imun,6) = temp_var(9)
!!$                    sav_inc(i3,imun) =  temp_var(12)
                 END IF

!!$                 IF (myrank_t == 0) write(*,'(a,5f20.10)') 'av_sav_2t, av_ab_2t', av_sav_2t, av_ab_2t, temp_var(15), temp_var(12), temp_var(14)                 
              END IF
              
           END DO
        END DO
     ENDDO
     av_sav_2t = av_sav_2t/DBLE(count_av)
     av_log_ab_2t = av_log_ab_2t/DBLE(count_av)
     std_log_ab_2t = (second_mom_log_ab_2t/DBLE(count_av) - av_log_ab_2t**2d0)**0.5d0
     DO ii = 1, Retage - 1
        prob_age_2t(ii) = DBLE(count_age(ii))/DBLE(count_av)
     END DO
     prob_edu_2t(1) = DBLE(count_low_edu)/DBLE(count_av)
     prob_edu_2t(2) = 1d0 - prob_edu_2t(1)
     prob_mun_2t(1) = DBLE(i1)/DBLE(i1+i2+i3)
     prob_mun_2t(2) = DBLE(i2)/DBLE(i1+i2+i3)
     prob_mun_2t(3) = DBLE(i3)/DBLE(i1+i2+i3)
          
     CLOSE(18)

!!$     ! Average savings of second term mayors in the base case
!!$     av_sav_2t = SUM(var_base(:,:,2,2), var_base(:,:,2,1) > 0.5d0)/DBLE(Nsimt*Nactmun)
!!$     av_ab_2t = SUM(var_base(:,:,2,3), var_base(:,:,2,1) > 0.5d0)/DBLE(Nsimt*Nactmun)

     IF (myrank_t == 0) write(*,'(a,3f20.10)') 'av_sav_2t, av_log_ab_2t, std_log_ab_2t', av_sav_2t, av_log_ab_2t, std_log_ab_2t
     IF (myrank_t == 0) write(*,'(a,7f20.10)') 'prob_age_2t', prob_age_2t
     IF (myrank_t == 0) write(*,'(a,2f20.10)') 'prob_edu_2t', prob_edu_2t
     IF (myrank_t == 0) write(*,'(a,3f20.10)') 'prob_mun_2t', prob_mun_2t

  END IF

  counter = 0
  DO j = 1, Nsimt

     counter = counter + 1
     IF ((counter <= ntasks*nprocs_t .AND. counter >= myrank_t*ntasks + 1 .AND. counter <= myrank_t*ntasks + ntasks) .OR. &
          (counter > ntasks*nprocs_t .AND. counter - ntasks*nprocs_t - 1 == myrank_t)) THEN
        
        DO i = 1, Nactmun

           CALL cpu_time(ts3)
           
           DO tt = 1, Totterms

              infeasible = 0
              
              thief = 0
              thief3 = 0
              aud3 = 0 
              auditf3 = 0

              uud(1:6) = udraw(i,j,tt,:)
              uud_age(1:2) = age_udraw(i,j,tt,:)
              
              !We multiply the zzpr_shocks distributed N(0,1) times the std deviation of the distribution to get a N(0,sigma2)
              zzpr_shock_t = sigma_zzpr*zzpr_shock(i,j,tt)
!!$              zzpr_shock_t = 0d0
              !Shock to the decision to run
              run_shock_t = sigma_run*rnormprob(i,j,tt) !zzpr_shock(i,j,tt)
!!$              if (myrank_t == 0) write(*,*) 'run_shock_t', run_shock_t, sigma_run, rnormprob(i,j,tt)

              !We transform the draw from the standard normal in draws from a log-normal(mu_fine,sigma_fine)
              !Fraction of the fine applied in case the mayor is audited this term
              frac_fine_p = frac_fine_c
              !Fraction of the fine used in mayordecisions to determine the amount of wealth he will have in the next term if he steals
              !frac_fine_p must always be equal frac_fine_n in the previous term: what the mayor believes will be applied if audited, must be equal to what will actually be applied if audited 

              ! We draw a probability to see if they mayor will be convicted in case s/he is audited and had stolen
              prob_conv_vec(:) = (/ prob_conv, 1d0-prob_conv /)                    
              CALL rmultinomial(uud(2), 2, prob_conv_vec, outfine)

              ! outfine = 1 if convicted, returns the funds stolen, and pays the fine
              ! outfine = 2 if not convicted, doesn't return the funds stolen, and doesn't pay the fine
              IF (outfine == 1) THEN
                 frac_fine_c = EXP(sigma_fine*fine_norm(i,j,tt) + mu_fine) 
              ELSE
                 frac_fine_c = -1d0
              END IF
!!$              IF (frac_fine_c > max_frac_fine) frac_fine_c = max_frac_fine

!!$              IF (flag_no_selec_sav == 1) THEN
!!$                 ! electoral outcome in the base case at the start of period tt
!!$                 outelect_base = var_base(j,i,tt,1)
!!$                 ! savings in the base case at the beginning of the period
!!$                 msaving_base = var_base(j,i,tt,2)
!!$              END IF

              IF (tt == 1) THEN

                 !In the first term set the prob of being audited to 0.04 following the data in Brazil
                 ipau = 1
                 ipauf = 1
                 ! In tt = 1, all mayors are term 1, so have no past
                 ipastf = 1
                 ! Since they have no past, their past stealing is equalt to zero
                 stealf = 0d0
                 ! And the previous incumbents decided not to run
                 d_run = 0d0

                 !Observed/Conditioning Variables
                 !imun descibes the type of municipality in terms of population; we have six types.
                 imun = INT(initdata(i,10))
                 subqcons = subqcons_vec(imun)
                 popsizeg = popsize(imun)
                 auditf = INT(initdata(i,2))
                 
                 mayorage = INT(initdata(i,3))
                 IF (flag_change_age_policy == 1) THEN
!!$                    IF (myrank == 0) write(*,'(a,1i5)') 'mayor age b', mayorage
                    pr_ch_age_vec(:) = (/ pr_ch_age, 1d0-pr_ch_age /)                    
                    CALL rmultinomial(uud_age(1), 2, pr_ch_age_vec, out_ch_age)
                    IF (out_ch_age == 1) mayorage = mayorage + ch_m_age
!!$                    IF (myrank == 0) write(*,'(a,3i5,3f20.10)') 'mayor age a', mayorage, ch_m_age, out_ch_age, uud_age(1), pr_ch_age_vec(:)
                 END IF
                 
                 ! For now, censoring the age of mayors: later we clean up the sample
                 mayorage_elec = MIN(Retage-1,mayorage)
                 incumb_age = mayorage_elec
                 ieducf = INT(initdata(i,4))

                 fundsf = per_funds*initdata(i,6)/norm
                 ! Since they have no past, their past fraction stolen is equalt to zero
                 fr_stolen = stealf/fundsf

!!$                 zzprf = initdata(i,5) * EXP(zzpr_shock_t)
                 zzprf = initdata(i,5)
                 nterms = INT(initdata(i,8))

                 !We don't need qconsf for Brasil since all mayors are in their first term in 1996
                 !However we need to give it a value for drawinitunob

                 !We compute EXP(qconsvec) to have the level
                 qconsf = qconsvec(imun,2)

                 !Wage of the mayor, it depends on the type of municipality
                 !We use the wage from our grid                 
!!$                 mwagef =  n_times_wage*initdata(i,7)/norm
                 mwagef =  wage_mayor(imun)

                 ! Potential wage in the private sector for the incumbent
                 potential_wage = wage_pastmayor(mayorage,imun,ieducf,ipastf)

                 !Outcome of the election in the first period, 1 means that the incumbent wins the election

                 !In the data, at tt = 1, all mayors are new mayors
                 outelect = INT(initdata(i,9))       
                 yy = 1d0

                 IF (outelect == 0 .AND. auditf == 1) WRITE(*,*) 'wrong initial data, outelect, auditf',  outelect, auditf

                 !We draw ability from a log-normal distribution
!!$                 ability = EXP(sigma_a*ability_norm(i,1,tt)+mu_a)
                 ability = EXP(sigma_a*ability_norm(i,j,tt)+mu_a)
!!$                 ability = sigma_a*ability_norm(i,j,tt)+mu_a

                 !If we fix ability at some level
                 IF (flag_fixed_ability == 1) ability = fixed_ability

                 elap = 1
                 IF (uud(6) >= pappeal(1)) THEN
                    elap = 2
                 END IF
                 
                 msavingf = initdata(i,18)/norm

                 !In the first period this is irrelevant because all mayors are in their 1st term
                 thief=0
!!$                 IF(stealf > l_steal_prob/norm) thief =1
                 IF(stealf/fundsf > l_fr_stolen) thief =1

                 !We save the variables for the current period
                 simout(i,tt,1) = DBLE(i) 
                 simout(i,tt,2) = DBLE(j) 
                 simout(i,tt,3) = DBLE(tt) 
                 simout(i,tt,4) = fr_stolen
                 !qconsf is in levels
                 simout(i,tt,5) = qconsf
                 simout(i,tt,6) = mwagef
                 simout(i,tt,7) = DBLE(auditf)
                 simout(i,tt,8) = stealf
                 simout(i,tt,9) = fundsf
!!$                 simout(i,tt,10) = zzprf/EXP(zzpr_shock_t)
                 simout(i,tt,10) = zzprf
                 simout(i,tt,11) = ability
                 
                 !This is only useful is we write out data
                 IF (flag_simulate_data == 1 .OR. (flag_policy_evaluation == 1 .AND. nboot==1)) THEN
                    simoutt(i,j,tt,1) = DBLE(i) 
                    simoutt(i,j,tt,2) = DBLE(j) 
                    simoutt(i,j,tt,3) = DBLE(tt) 
                    simoutt(i,j,tt,4) = fr_stolen
                    simoutt(i,j,tt,5) = qconsf
                    simoutt(i,j,tt,6) = mwagef
                    simoutt(i,j,tt,7) = DBLE(auditf)
                    simoutt(i,j,tt,8) = stealf
                    simoutt(i,j,tt,9) = fundsf
!!$                    simoutt(i,j,tt,10) = zzprf/EXP(zzpr_shock_t)
                    simoutt(i,j,tt,10) = zzprf
                    simoutt(i,j,tt,11) = ability
                 END IF

              ELSE !IF tt==1

                 IF (flag_no_selec_age == 1 .OR. flag_no_selec_educ == 1 .OR. flag_no_selec_ab == 1 .OR. flag_no_selec_sav == 1 .OR. flag_no_selec_funds == 1 .OR. flag_no_selec_zzpr == 1 .OR. flag_no_selec_mun == 1) THEN

                    CALL RANDOM_NUMBER(uselect)
                    IF (flag_no_selec_mun == 1) THEN
                       CALL rmultinomial(uud(1), Nmun, prob_mun_2t, outmun)
                       IF (outmun == 1) THEN
                          iselect = 1 + FLOOR(i1*uselect) 
                       ELSE IF (outmun == 2) THEN
                          iselect = 1 + FLOOR(i2*uselect) 
                       ELSE IF (outmun == 3) THEN
                          iselect = 1 + FLOOR(i3*uselect) 
                       END IF
                       imun = outmun
                       popsizeg = popsize(imun)
                       subqcons = subqcons_vec(imun)
                       mwagef =  wage_mayor(imun)
                    ELSE
                       IF (imun == 1) THEN
                          iselect = 1 + FLOOR(i1*uselect) 
                       ELSE IF (imun == 2) THEN
                          iselect = 1 + FLOOR(i2*uselect) 
                       ELSE IF (imun == 3) THEN
                          iselect = 1 + FLOOR(i3*uselect) 
                       END IF
                    END IF

!!$                    if (myrank_t == 0) write(*,'(a,3i4,3f20.10)') 'random draw', outnum, imun, iselect, popsizeg, subqcons, mwagef
                    
                 END IF
                 
                 !In the first term set the prob of being audited to 0.04 following the data in Brazil
                 IF (tt == 2) THEN
                    !ipau is used to draw the audit outcome at the beginning of the term for the mayor in power the previous term this is why we have ipau = 1 fro tt = 2
                    ipau = 1
                    !ipauf is used for the betas to construct value functions for future periods
                    ipauf = 2
                 ELSE
                    ipau = 2
                    ipauf = 2                    
                 END IF
                 
                 !If nterms >= Nnterms the mayor must step down
                 dft = 0
                 IF (nterms >= Nnterms) dft = 1
 
                 !Incumbent starting variables
                 stealf = newstealf
                 ! Both stealf and funds are for the previous period
                 fr_stolen = stealf/fundsf
                 !The following index is needed to choose the betas for a mayor who stole more than x%
                 index_steal_x = 1
                 IF (fr_stolen > l_fr_stolen_x) index_steal_x = 2                 
                 !qconsf is in levels
                 qconsf = newqconsf
                 msavingf = mnewsavingf                 

                 ieducf = ieducf
                 incumb_educ = ieducf
                 mayorage = mayorage_elec + 1
                 incumb_age = mayorage

                 !audit
                 pauditt(1:2) = paudit(0:1,ipau)
!!$                 IF (flag_audit_nterm == 1 .AND. (nterms == Nnterms .OR. mayorage >= Retage)) pauditt(1:2) = paudit_nterm(0:1)
                 IF (flag_audit_nterm == 1 .AND. nterms == Nnterms) pauditt(1:2) = paudit_nterm(0:1)
                 CALL rmultinomial(uud(5), 2, pauditt, aud)
                 auditf = aud-1

!!$                 if (myrank == 0 .AND. tt == 3) write(*,'(a,1i4,2f20.10)') 'pauditt(1:2)', auditf, pauditt(1:2)
                 
                 !Policy that audits all mayors that were audited and caught stealing in the past
                 !This is the right place for setting auditf = 1, beacause the policy has an effect only if the mayor was caught stealing in the past and not in the current period
                 !Notice that the policy has an effect only through the value function if Nnterms<=2, because only mayors at the end of the second term can have ipastf = 3
!!$                 !If Nnterms>2, then it can also have a direct effect
!!$                 IF (flag_policy_audit == 1) THEN
!!$                    IF (stealf > 0d0) auditf = 1
!!$                 END IF

                 ! We draw a challenger facing the incumbent in the election or becoming the incumbent if the current mayor does not run
                 ! We have to draw the challenger here and not later in case the incumbent must step down because nterms > Nnterms
                 ! The challenger has three observable characteristics: age, wealth, and education                 

                 ! When we draw the challenger we should draw from nterms = 1
                 IF (imun == 1) THEN

                    aenoc_mun = INT(uud(4) * (DBLE(countNnterms_mun1(1))-1.0d0)) + 1                       
                    ch_age = INT(nonpardistnoc_mun1(aenoc_mun,1,1))
                    ch_educ = INT(nonpardistnoc_mun1(aenoc_mun,1,2))
                    ch_saving = nonpardistnoc_mun1(aenoc_mun,1,5)
                    
                 ELSE IF (imun == 2) THEN
                    
                    aenoc_mun = INT(uud(4) * (DBLE(countNnterms_mun2(1))-1.0d0)) + 1
                    ch_age = INT(nonpardistnoc_mun2(aenoc_mun,1,1))
                    ch_educ = INT(nonpardistnoc_mun2(aenoc_mun,1,2))
                    ch_saving = nonpardistnoc_mun2(aenoc_mun,1,5)
                    
                 ELSE IF (imun == 3) THEN
                    
                    aenoc_mun = INT(uud(4) * (DBLE(countNnterms_mun3(1))-1.0d0)) + 1                       
                    ch_age = INT(nonpardistnoc_mun3(aenoc_mun,1,1))
                    ch_educ = INT(nonpardistnoc_mun3(aenoc_mun,1,2))
                    ch_saving = nonpardistnoc_mun3(aenoc_mun,1,5)
                    
                 END IF

                 IF (flag_change_age_policy == 1) THEN
!!$                    IF (myrank == 0) write(*,'(a,1i5)') 'challenger age b', ch_age
                    pr_ch_age_vec(:) = (/ pr_ch_age, 1d0-pr_ch_age /)                    
                    CALL rmultinomial(uud_age(2), 2, pr_ch_age_vec, out_ch_age)
                    IF (out_ch_age == 1) ch_age = ch_age + ch_m_age
!!$                    IF (myrank == 0) write(*,'(a,3i5,3f20.10)') 'challenger age a', ch_age, ch_m_age, out_ch_age, uud_age(2), pr_ch_age_vec(:)
                 END IF
                 
!!$                 !For the second term we observe in the data, we use mayor's age and education from the data
!!$                 IF (tt == 2) THEN
!!$                    
!!$                    ch_age = INT(initdata(i,15))
!!$                    ch_educ = INT(initdata(i,16))
!!$                    
!!$                 END IF
                 
                 IF (flag_no_selec_age == 1) THEN
                    ch_age = INT(var_base(iselect,imun,3))
                 END IF
                 IF (flag_no_selec_educ == 1) THEN
                    ch_educ = INT(var_base(iselect,imun,4))
                 END IF
                 
                 IF (flag_no_selec_sav == 1) THEN
                    ch_saving = var_base(iselect,imun,1)
                 END IF
                 
                 ch_age = MIN(Retage-1,ch_age)
                 
                 !If the mayor has served the maximum number of terms or s/he is older than the maximum age (Retage)
                 !the challenger automatically wins the election.
                 IF (dft == 1 .OR. mayorage >= Retage .OR. infeasible == 1) THEN
                    yy = -99999999999d0
                    outelect = -999999
                    infeasible = 0
                    d_run = -1d0
                    GO TO 80
                 END IF

                 IF (infeasible == 1) THEN
                    yy = -77777777777777d0
                    outelect = -777777
                    infeasible = 0
                    d_run = -777d0
                    GO TO 80
                 END IF

                 !ln_relqcons is first used in the computation of the value function
                 !the value function depends on log of public consumption acounting for the degree of rivalry: qconsvec = log(Q/P**eta)
!!$                 ln_relqcons = LOG(qconsf)-eta*LOG(popsizeg)
                 ln_relqcons = LOG(qconsf)-LOG(popsizeg)
                 
                 !The first decision of the incumbent is whether to run or not.
                 !We first compute the mayor's value function if s/he decides to run this period 
                 !ALL THIS IS CONDITIONAL ON KNOWING THE RESULT OF AUDIT!
                 !For now, this is only valid with Nnterms == 2; it has to be modified to allow for iterm =3 or more!

!!$                 IF (msavingf < 0d0) write(*,'(a,1f20.10)') 'neg sav 1 1', msavingf

                 ! In the following, ipastf = 1 is no audit, = 2 if audit
                 msaving_bf = msavingf
                 IF(auditf == 0) THEN
                    !In this case, we keep the ipastf from the previous term to use the correct wage for the tax policy wage
                    rundec_netW = msavingf
                    ipastf = 1
                 ! fr_stolen is for the previous period
                 ELSE IF(auditf == 1 .AND. fr_stolen <= l_fr_stolen) THEN
                    ipastf = 2
                    rundec_netW = msavingf
                 ! fr_stolen is for the previous period
                 ELSE IF(auditf == 1 .AND. fr_stolen > l_fr_stolen) THEN
                    !We have to use the frac_fine saved from the previous term
                    frac_fine_g = frac_fine_p
                    rundec_netW = msavingf - fine_fun(stealf)
                    msavingf = rundec_netW
                    ipastf = 2
                 END IF

                 d_audit = 0d0
                 IF (auditf == 1) d_audit = 1d0
                 
                 ! Potential wage in the private sector for the incumbent
                 ! For past mayors we always have ipastf = 1 (the last variables in wage_pastmayor), ipastf = 2 is for voters
                 potential_wage = wage_pastmayor(mayorage,imun,ieducf,1)

!!$                 will_to_pay = 0d0
!!$                 ! WE COMPUTE THE WILLINGNESS TO PAY HERE USING NTERMS + 1, IT IS THE VALUE BEFORE THE DECISION TO RUN AND THE ELECTIONS
!!$                 ! We compute the value of the policy
!!$                 IF (flag_simulate_data == 1 .AND. flag_base_policy == 0) THEN
!!$                                        
!!$                    ! We should compute the value function for the median resident only in the first period we consider
!!$                    ! Since we compute the value function, it should include the effect of the policy for the rest of the median resident's life
!!$                    ! THIS IS NOT TRUE BECAUSE IN tt == 1 ALL MAYORS ARE FIRST TERM AND THEY MAKE DIFFERNET DECISIONS
!!$                    ! We use median age, education, and wealth conditional on the size of the municipality
!!$                    med_edu_i = med_edu(imun)
!!$                    med_age_g = med_age(imun)
!!$                    med_ipast_i = med_ipast(imun)
!!$                    med_wealth_g = med_wealth(imun)
!!$                    pmwage_g = wage_vot(imun)
!!$                    
!!$                    IF (ABS(delta-1d0) < small_no) THEN
!!$                       ELNqcons = LOG(qconsf/((popsizeg)**eta))
!!$                    ELSE
!!$                       ELNqcons = (qconsf/((popsizeg)**eta))**(1d0-delta)/(1d0-delta)
!!$                    END IF
!!$                    ELNqcons_g =LOG(qconsf/((popsizeg)**eta))
!!$                 
!!$                    !Find the point on the saving domain where we need to compute the approximated value function
!!$                    DO isav = 1, Nwealth
!!$                       IF (msavingf <  wealth(isav)) THEN
!!$                          ub_savi = isav
!!$                          EXIT
!!$                       END IF
!!$                       ub_savi = Nwealth
!!$                    END DO
!!$                    
!!$                    ALLOCATE(betat(MMage_steal_q),betat_base(MMage_steal_q))
!!$                    betat = 0d0
!!$                    betat_base = 0d0
!!$                    MMg = Nwealth
!!$                    elapv = 1
!!$                    IF (pappeal(1) < 0.5d0) elapv = 2
!!$                    
!!$                    ! VOTER'S VALUE FOR THE BASE CASE
!!$                    betat_base(:) = betapmayor_base(med_age_g,imun,nterms+1,ipastf,ieducf,med_ipast_i,ub_savi,ipauf,med_edu_i,elapv,:)
!!$                    
!!$                    IF (ub_sav_med == 1) THEN
!!$                       
!!$                       multi_slope = DBLE(CEILING(ABS((med_wealth_g - wealth(1)))/inc_wealth))
!!$                       V_slope_base = (1d0+multi_slope)*(betat_base(2) - betat_base(1))/(wealth(2) - wealth(1))
!!$                       
!!$                    ELSE
!!$                       
!!$                       V_slope_base = (betat_base(ub_sav_med) - betat_base(ub_sav_med-1))/(wealth(ub_sav_med) - wealth(ub_sav_med-1))
!!$                       
!!$                    END IF
!!$                    
!!$                    const_base = betat_base(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+1) + &
!!$                         betat_base(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+2)*incumb_age + &
!!$                         betat_base(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+3)*incumb_age**2d0 + &
!!$                         betat_base(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+4)*fr_stolen + &
!!$                         betat_base(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+5)*ELNqcons_g + betat_base(ub_sav_med)
!!$                    
!!$                    value_base = const_base + V_slope_base*(med_wealth_g - wealth(ub_sav_med))
!!$                    
!!$                    ! VOTER'S VALUE FOR THE POLICY CASE
!!$                    betat(:) = betapmayor(med_age_g,imun,nterms+1,ipastf,ieducf,med_ipast_i,ub_savi,ipauf,med_edu_i,elapv,:)
!!$                    IF (ub_sav_med == 1) THEN
!!$                       
!!$                       multi_slope = DBLE(CEILING(ABS((med_wealth_g - wealth(1)))/inc_wealth))
!!$                       V_slope = (1d0+multi_slope)*(betat(2) - betat(1))/(wealth(2) - wealth(1))
!!$                       
!!$                    ELSE
!!$                       
!!$                       V_slope = (betat(ub_sav_med) - betat(ub_sav_med-1))/(wealth(ub_sav_med) - wealth(ub_sav_med-1))
!!$                       
!!$                    END IF
!!$                    
!!$                    const_voter = betat(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+1) + &
!!$                         betat(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+2)*mayorage_elec + &
!!$                         betat(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+3)*mayorage_elec**2d0 + &
!!$                         betat(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+4)*fr_stolen + &
!!$                         betat(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+5)*ELNqcons_g
!!$                    
!!$                    value_voter = const_voter + V_slope*(med_wealth_g - wealth(ub_sav_med)) + betat(ub_sav_med)
!!$                 
!!$                    ! WE FIND THE LEVEL OF WEALTH THAT MAKES THE INDIVIDUAL INDIFFERENT BETWEEN THE POLICY AND BASE CASE
!!$                    ! WE HAVE TO USE THE BETAS FOR THE POLICY
!!$                    DO isav = 1, Nwealth
!!$                       IF (const_voter + betat(isav) > value_base) THEN
!!$                          ub_sav_indiff = isav
!!$                          EXIT
!!$                       END IF
!!$                       ub_sav_indiff = Nwealth
!!$                    END DO
!!$                    
!!$                    IF (ub_sav_indiff == 1) THEN
!!$                       
!!$                       multi_slope = DBLE(CEILING(ABS((med_wealth_g - wealth(1)))/inc_wealth))
!!$                       V_slope_indiff = (1d0+multi_slope)*(betat(2) - betat(1))/(wealth(2) - wealth(1))
!!$                       
!!$                    ELSE
!!$                       
!!$                       V_slope_indiff = (betat(ub_sav_indiff) - betat(ub_sav_indiff-1))/(wealth(ub_sav_indiff) - wealth(ub_sav_indiff-1))
!!$                       
!!$                    END IF
!!$                    
!!$                    const_indiff = const_voter + betat(ub_sav_indiff)
!!$                    
!!$                    !write(*,'(a,15f20.10)') 'beta_base 0', betat_base(:)
!!$                    !write(*,'(a,15f20.10)') 'beta_base 1', betapmayor_base(med_age_g,imun,nterms+1,ipastf,ieducf,med_ipast_i,ub_savi,ipauf,med_edu_i,elapv,:)
!!$                    !write(*,'(a,15f20.10)') 'betat', betat(:)
!!$                    !write(*,'(a,15f20.10)') 'betat', betapmayor(med_age_g,imun,nterms+1,ipastf,ieducf,med_ipast_i,ub_savi,ipauf,med_edu_i,elapv,:)
!!$                    
!!$                    !write(*,'(a,2i4,9f20.10)') 'will_to_pay, med_wealth_g, wealth_indiff', ub_sav_med, ub_sav_indiff, V_slope_base, V_slope, V_slope_indiff, med_wealth_g, value_voter, value_base
!!$                    !write(*,'(a,11i5)') 'indeces', med_age_g,imun,nterms+1,ipastf,ieducf,med_ipast_i,ub_savi,ipauf,med_edu_i,elapv
!!$                    
!!$                    wealth_indiff = (const_base - const_indiff + V_slope_base*(med_wealth_g - wealth(ub_sav_med)) + V_slope_indiff*wealth(ub_sav_indiff))/V_slope_indiff
!!$                    value_indiff = const_indiff + V_slope_indiff*(wealth_indiff - wealth(ub_sav_indiff))
!!$                    
!!$                    IF (ABS(value_base - value_indiff) >= 0.0001d0) write(*,'(a,2i4,8f14.7)') 'wrong calculations wtp', ub_sav_med, ub_sav_indiff, value_base, V_slope_indiff*(wealth_indiff - wealth(ub_sav_indiff)) + betat(ub_sav_indiff), value_voter, med_wealth_g, wealth_indiff, betat_base(ub_sav_med), betat(ub_sav_med), betat(ub_sav_indiff)
!!$                    
!!$                    will_to_pay = med_wealth_g - wealth_indiff
!!$
!!$                    !if (myrank_t == 0) THEN
!!$                       !write(*,'(a,2i4,9f20.10)') 'will_to_pay, med_wealth_g, wealth_indiff', ub_sav_med, ub_sav_indiff, V_slope_base, V_slope, V_slope_indiff, will_to_pay, med_wealth_g, wealth_indiff, value_voter, value_base, value_indiff
!!$                       !write(*,'(a,15f20.10)') 'beta_base', betat_base(:)
!!$                       !write(*,'(a,15f20.10)') 'betat', betat(:) 
!!$                    !END if
!!$                    
!!$                    !if (myrank_t == 0) write(*,'(a,11i5)') 'indeces', med_age_g,imun,nterms+1,ipastf,ieducf,med_ipast_i,ub_savi,ipauf,med_edu_i,elapv
!!$
!!$                    av_will_to_pay_t(i,tt) = av_will_to_pay_t(i,tt) + will_to_pay                    
!!$                    
!!$                    DEALLOCATE(betat,betat_base)
!!$                    
!!$                 END IF

                 ALLOCATE(betart(MMability_steal_q))
                 betart =  betasimrun(mayorage,imun,ipastf,nterms+1,ipauf,ieducf,elap,:)
                 
                 !Find the point on the saving domain where we need to compute the approximated value function
                 DO isav = 1, Nwealth
                    IF (rundec_netW <  wealth(isav)) THEN
                       ub_sav = isav
                       EXIT
                    END IF
                    ub_sav = Nwealth
                 END DO
                       
                 !COMPUTE THE LINEAR APPROXIMATION
                 IF (ub_sav == 1) THEN
                          
                    multi_slope = DBLE(CEILING(ABS((rundec_netW - wealth(1)))/inc_wealth))
                    V_slope = (1d0+multi_slope)*(betart(2) - betart(1))/(wealth(2) - wealth(1))
                    
                 ELSE
                    
                    V_slope = (betart(ub_sav) - betart(ub_sav-1))/(wealth(ub_sav) - wealth(ub_sav-1))
                    
                 END IF

                 ! We need to use the frac_fine that affects stealing in the previoius term: frac_fine_p
                 emaxr = betart(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+1) + betart(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+2)*ability + &
                      betart(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+3)*ability**2d0 + betart(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+4)*fr_stolen*d_audit + &
                      betart(MMability_steal_q-Npol_ab-Npol_st-Npol_q-1+5)*ln_relqcons + V_slope*(rundec_netW - wealth(ub_sav)) + betart(ub_sav)

                 emaxr = emaxr + run_shock_t - cost_run_frac
                 
!!$                 ! If we evaluate the no-run policy it is easier to analyse the data without the shock
!!$                 IF (flag_no_run == 1 .AND. auditf == 1 .AND. fr_stolen > l_fr_stolen) emaxr = emaxr - run_shock_t
!!$                 IF (flag_no_run_x == 1 .AND. auditf == 1 .AND. fr_stolen > l_fr_stolen_x) emaxr = emaxr - run_shock_t

                 !We then compute the mayor's value function if s/he decides not to run this period 
                 !This is the expected value for this period. No exogenous variables (other than audit) is realized
                 ALLOCATE(betanort(Nwealth))
                 betanort = betanornew(mayorage,imun,ipastf,ipauf,ieducf,elap,:)

                 !COMPUTE THE LINEAR APPROXIMATION
                 IF (ub_sav == 1) THEN
                    
                    multi_slope = DBLE(CEILING(ABS((rundec_netW - wealth(1)))/inc_wealth))
                    V_slope = (1d0+multi_slope)*(betanort(2) - betanort(1))/(wealth(2) - wealth(1))
                    
                 ELSE
                    
                    V_slope = (betanort(ub_sav) - betanort(ub_sav-1))/(wealth(ub_sav) - wealth(ub_sav-1))
                    
                 END IF
                 
                 emaxnor = V_slope*(rundec_netW - wealth(ub_sav)) + betanort(ub_sav)

                 dum_steal = 0d0
                 dum_steal_x = 0d0
                 IF (fr_stolen > l_fr_stolen) dum_steal = 1d0
                 IF (fr_stolen > l_fr_stolen_x) dum_steal_x = 1d0
                 
                 IF (emaxnor > emaxr) THEN
                    yy = -111111111d0
                    outelect = -111111
                    !Dummy to record whether the incumbent run for reelection
                    d_run = 0d0
                    
                 ! If we want to implement a policy in which the mayor looses the election with probability 1 if he steals (equivalent to not allowing the mayor to run for reelection)
                 ELSE IF (flag_no_run == 1 .AND. DBLE(auditf)*dum_steal > 0.5d0) THEN
                    d_run = 0d0
                    yy = -1000d0
                    outelect = 0
                 ! If we want to implement a policy in which the mayor looses the election with probability 1 if he steals more than x (equivalent to not allowing the mayor to run for reelection)
                 ELSE IF (flag_no_run_x == 1 .AND. DBLE(auditf)*dum_steal_x > 0.5d0) THEN
                    d_run = 0d0
                    yy = -1000d0
                    outelect = 0
                 ELSE IF (flag_norun_prob_policy == 1 .AND. DBLE(auditf)*dum_steal_x > 0.5d0) THEN
                    d_run = 0d0
                    yy = -1000d0
                    outelect = 0                    
                 ELSE
                    
                    !Dummy to record whether the incumbent run for reelection
                    d_run = 1d0
                                           
                    ALLOCATE(beta_i(MMage_steal_q), beta_c(MMage))
                    beta_i = 0d0
                    beta_c = 0d0
                    
                    ! In the data, we only have ieduv = 1
                    ieduv = 1

                    !For the challenger, find the point on the saving domain where we need to compute the approximated vemax_voter_c
                    DO isav = 1, Nwealth
                       IF (ch_saving < wealth(isav)) THEN
                          ub_sav_ch = isav
                          EXIT
                       END IF
                       ub_sav_ch = Nwealth
                    END DO
                    
                    ! We normalize popsize to 1
                    popsizet = 1d0
                    ! Dummy for being audited
                    d_past = 0d0
                    IF (ipastf == 2) d_past = 1d0
                    alpha_nstt = 0d0
                    IF (ipastf == 2 .AND. fr_stolen <= 0d0) alpha_nstt = alpha_nst
                    alpha_eat = 0d0
                    IF (ipastf == 2 .AND. fr_stolen > 0d0) alpha_eat = - alpha_ea
!!$                    IF (ipastf == 2 .AND. fr_stolen > 0d0) alpha_eat = rnorm(i,j,tt) - alpha_ea
!!$                    IF (ipastf == 2 .AND. fr_stolen > 0d0) alpha_eat = alpha_ea*rnorm(i,j,tt)              
!!$                    IF (elap == 2) alpha_eat = alpha_ea
                    av_qcons = SUM(av_qcons_vec(1:Retage-1,imun,ipauf))/DBLE(Retage-1)
!!$                    diff_qcons = ln_relqcons - LOG(av_qcons)
!!$                    diff_qcons = ln_relqcons
!!$                    ! qconsvec(imun,2) is average per capita public consumption in municipality imun in the data
                    diff_qcons = ln_relqcons - LOG(qconsvec(imun,2))
                    voter_age = 5
                          
                    ! We only know the reelection probability for discrete value of ability and wealth
                    ! We therefore use the approximations of the value function for voters if an incumbent or a challenger is elected computed in emaxfun.f90
                    ! We use the voter of age 5 as the average voter
                    beta_i(:) = beta_voter_i(voter_age,imun,nterms+1,ipastf,ieducf,ipauf,ieduv,elap,:)

                    ! We construct the approximated value function if the incumbent wins
                    IF (ub_sav == 1) THEN
                       
                       multi_slope = DBLE(CEILING(ABS((rundec_netW - wealth(1)))/inc_wealth))
                       V_slope = (1d0+multi_slope)*(beta_i(2) - beta_i(1))/(wealth(2) - wealth(1))
                       
                    ELSE
                       
                       V_slope = (beta_i(ub_sav) - beta_i(ub_sav-1))/(wealth(ub_sav) - wealth(ub_sav-1))
                       
                    END IF

                    emax_voter_i = beta_i(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+1) + beta_i(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+2)*DBLE(mayorage) + &
                         beta_i(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+3)*DBLE(mayorage)**2d0 + beta_i(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+4)*fr_stolen*d_audit + &
                         beta_i(MMage_steal_q-Npol_age-Npol_st-Npol_q-1+5)*ln_relqcons + V_slope*(rundec_netW - wealth(ub_sav)) + beta_i(ub_sav)

                    ! We construct the approximated value function if the incumbent wins
                    beta_c(:) = beta_voter_c(voter_age,imun,ch_educ,ipauf,ieduv,elap,:)

                    IF (ub_sav_ch == 1) THEN
                       
                       multi_slope = DBLE(CEILING(ABS((ch_saving - wealth(1)))/inc_wealth))
                       V_slope = (1d0+multi_slope)*(beta_c(2) - beta_c(1))/(wealth(2) - wealth(1))
                       
                    ELSE
                       
                       V_slope = (beta_c(ub_sav_ch) - beta_c(ub_sav_ch-1))/(wealth(ub_sav_ch) - wealth(ub_sav_ch-1))
                       
                    END IF

                    emax_voter_c = beta_c(MMage-Npol_age-1+1) + beta_c(MMage-Npol_age-1+2)*DBLE(ch_age) + beta_c(MMage-Npol_age-1+3)*DBLE(ch_age)**2d0 + &
                         V_slope*(ch_saving - wealth(ub_sav_ch)) + beta_c(ub_sav_ch)

                    arg_exp = -(emax_voter_i - emax_voter_c + alpha_ie + alpha_nstt + alpha_eat + alpha_st*d_past*fr_stolen + alpha_q*diff_qcons)

                    ! We compute the probability that a voter with education ieduv votes for the incumbent
                    IF (arg_exp > -50d0 .AND. arg_exp < 50d0) THEN
                       pp = 1d0/(1d0 + EXP(arg_exp))
                    ELSE IF (arg_exp <= -50d0) THEN
                       pp = 1d0-small_no
                    ELSE IF (arg_exp > 50d0) THEN
                       pp = small_no
                    END IF

                    Phi_arg = popsizet**0.5d0*(0.5d0 -pp)/(pp*(1d0-pp))**0.5d0
                       
                    Pr_win_i = 1d0 - pnorm(Phi_arg)

                    Pwin = (/ 1d0 - Pr_win_i, Pr_win_i /)
                    
                    CALL rmultinomial(uud(3), 2, Pwin, outelect)
                    outelect = outelect - 1
                    yy = 1d0
                    IF (outelect == 0) yy = 0d0

!!$                    IF (myrank == 0) write(*,'(a,12f20.10)') 'Pwin', Pwin, yy, emax_voter_i, emax_voter_c, alpha_ie, alpha_nstt, alpha_eat, alpha_st*d_past*fr_stolen, alpha_q*diff_qcons
                    
                    DEALLOCATE(beta_i, beta_c)
                    
                 END IF
                 
                 DEALLOCATE(betart,betanort)
                 
                 !We save the variables for the last period
80               simout(i,tt,1) = DBLE(i)
                 simout(i,tt,2) = DBLE(j)
                 simout(i,tt,3) = DBLE(tt)
                 simout(i,tt,4) = fr_stolen
                 !newqconsf is in levels
                 simout(i,tt,5) = qconsf
                 simout(i,tt,6) = mwagef
                 simout(i,tt,7) = DBLE(auditf)
                 simout(i,tt,8) = stealf
                 !These two variables are saved later in the file
!!$                 simout(i,tt,9) = fundsf
!!$                 simout(i,tt,10) = zzprf/EXP(zzpr_shock_t)
                 simout(i,tt,11) = ability
                 
                 !This is only useful is we write out data
                 IF (flag_simulate_data == 1 .OR. (flag_policy_evaluation == 1 .AND. nboot==1)) THEN
                    simoutt(i,j,tt,1) = DBLE(i) 
                    simoutt(i,j,tt,2) = DBLE(j) 
                    simoutt(i,j,tt,3) = DBLE(tt) 
                    simoutt(i,j,tt,4) = fr_stolen
                    !newqconsf is in levels
                    simoutt(i,j,tt,5) = qconsf
                    simoutt(i,j,tt,6) = mwagef
                    simoutt(i,j,tt,7) = DBLE(auditf)
                    simoutt(i,j,tt,8) = stealf
                    !These two variables are saved later in the file
!!$                    simoutt(i,j,tt,9) = fundsf
!!$                    simoutt(i,j,tt,10) = zzprf/EXP(zzpr_shock_t)
                    simoutt(i,j,tt,11) = ability
                 END IF
                 
                 thief=0
                 
                 ! fr_stolen is for the previous period
                 IF(fr_stolen > l_fr_stolen) thief =1

                 ! This corresponds to what we have in probmayor: Phi(Xb) =  N_cdf(Xb <= epsilon); hence, we need yy = Xb - epsilon >= 0 (all cases in which -epsilon is small enough to be smaller than Xb
                 IF (outelect == 1) THEN
                    
                    !The incumbent stays in power
                    nterms = nterms + 1
                    mayorage_elec = mayorage

                    IF (flag_no_change_ability == 1) ability = EXP(sigma_a*ability_norm(i,j,tt)+mu_a)

!!$                    IF (myrank_t == 0) write(*,'(a,3i5,7f20.10,1f40.10)') 'ability in 2',i,j,tt, ability, sigma_a, ability_norm(i,j,tt), mu_a, d_run, emaxnor, emaxr, yy
                 
                 ELSE

                    !The challenger is elected
                    auditf = 0
                    ipastf = 1
                    nterms = 1

                    !The challenger's observable characteristics are drawn before the elections
                    mayorage = ch_age
                    ieducf = ch_educ
                    msavingf = ch_saving

                    !We draw ability from a log-normal distribution
                    ability = EXP(sigma_a*ability_norm(i,j,tt)+mu_a)

!!$                    IF (flag_no_selec_ab == 1) ability = EXP(std_log_ab_2t*ability_norm(i,j,tt)+av_log_ab_2t)
                    IF (flag_no_selec_ab == 1) ability = var_base(iselect,imun,2)
                    
                    !If we fix ability at some level
                    IF (flag_fixed_ability == 1) ability = fixed_ability
                    
!!$                    IF (myrank_t == 0) write(*,'(a,3i5,7f20.10,1f40.10)') 'ability cha 1',i,j,tt, ability, sigma_a, ability_norm(i,j,tt), mu_a, d_run, emaxnor, emaxr, yy

                    elap = 1
                    IF (uud(6) >= pappeal(1)) THEN
                       elap = 2
                    END IF
                    
                 END IF

                 !For the second term we observe in the data, we use funds and the private inputs from the data
                 !In tt == 2 a mayor can never be in the third term because in tt == 1 they all start as first term mayors
                 IF (tt == 2) THEN
                    
                    fundsf = per_funds*initdata(i,14)/norm
                    zzprf = initdata(i,17)
                    
                 ELSE 
                    
                    IF (imun == 1) THEN
                       
                       aenoc_mun = INT(uud(4) * (DBLE(Ntotmunnoc_mun1)-1.0d0)) + 1
                       !Funds provided by central governement
                       fundsf = per_funds*nonpardistnoc_uncond_mun1(aenoc_mun,1)
                       
                       !Private gdp
                       zzprf = per_funds*nonpardistnoc_uncond_mun1(aenoc_mun,2)
                       
!!$                       write(*,'(a,1i6,3f18.8)') 'funds', Ntotmunnoc_mun1, uud(4), nonpardistnoc_uncond_mun1(aenoc_mun,1), nonpardistnoc_uncond_mun1(aenoc_mun,2)
                       
                    ELSE IF (imun == 2) THEN
                       
                       aenoc_mun = INT(uud(4) * (DBLE(Ntotmunnoc_mun2)-1.0d0)) + 1
                       !Funds provided by central governement
                       fundsf = per_funds*nonpardistnoc_uncond_mun2(aenoc_mun,1)
                       
                       !Private gdp
                       zzprf = per_funds*nonpardistnoc_uncond_mun2(aenoc_mun,2)
                       
!!$                       write(*,'(a,1i6,3f18.8)') 'funds', Ntotmunnoc_mun2, uud(4), nonpardistnoc_uncond_mun2(aenoc_mun,1), nonpardistnoc_uncond_mun2(aenoc_mun,2)
                       
                    ELSE IF (imun == 3) THEN
                       
                       aenoc_mun = INT(uud(4) * (DBLE(Ntotmunnoc_mun3)-1.0d0)) + 1
                       !Funds provided by central governement
                       fundsf = per_funds*nonpardistnoc_uncond_mun3(aenoc_mun,1)
                       
                       !Private gdp
                       zzprf = per_funds*nonpardistnoc_uncond_mun3(aenoc_mun,2)
                       
!!$                       write(*,'(a,1i6,3f18.8)') 'funds', Ntotmunnoc_mun3, uud(4), nonpardistnoc_uncond_mun3(aenoc_mun,1), nonpardistnoc_uncond_mun3(aenoc_mun,2)
                       
                    END IF
                 END IF

                 IF (flag_no_selec_funds == 1) fundsf = var_base(iselect,imun,6)
                 IF (flag_no_selec_zzpr == 1) zzprf = var_base(iselect,imun,5)
                
                 !WE SAVE funds AND private inputs HERE BECAUSE THE VARIABLES ARE DRAWN FROM A DISTRIBUTION THAT DEPENDS ON nterms WHICH IS DETERMINED AFTER THE ELECTIONS
                 simout(i,tt,9) = fundsf
!!$                 simout(i,tt,10) = zzprf/EXP(zzpr_shock_t)
                 simout(i,tt,10) = zzprf
                 
                 !This is only useful if we write out data
                 IF (flag_simulate_data == 1 .OR. (flag_policy_evaluation == 1 .AND. nboot==1)) THEN
                    simoutt(i,j,tt,9) = fundsf
!!$                    simoutt(i,j,tt,10) = zzprf/EXP(zzpr_shock_t)
                    simoutt(i,j,tt,10) = zzprf
                 END IF
                 
              END IF  !END IF: Time dependent stuff!
              
              IF (mayorage_elec < Retage -1 .AND. nterms < Nnterms) THEN
                 ALLOCATE(betat1(MMability_steal_q),betat2(MMability_steal_q))
                 !inside mayorsolution.f90 we need the betas for t+1 for mayors with different types of stealing behavior the current period
                 IF (Npast > 1) THEN
                    betat1(:) = betanewcr(mayorage_elec+1,imun,1,nterms+1,ipauf,ieducf,elap,:)
                    betat2(:) = betanewcr(mayorage_elec+1,imun,2,nterms+1,ipauf,ieducf,elap,:)
                    IF (flag_no_run == 1 .OR. flag_norun_prob_policy == 1) THEN
                       ALLOCATE(betat3(Nwealth))
                       betat3(:) = betanornew(mayorage_elec+1,imun,2,ipauf,ieducf,elap,:)
                    END IF
                 ELSE
                    betat1(:) = betanewcr(mayorage_elec+1,imun,1,nterms+1,ipauf,ieducf,elap,:)
                 END IF
                 MMg = MMability_steal_q

              ELSE IF(mayorage_elec == Retage -1 .OR. nterms >=Nnterms) THEN 
                 ALLOCATE(betat1(Nwealth),betat2(Nwealth))
                 !If they cannot run, it is indifferent whether they steal in the current period
                 IF (Npast > 1) THEN
                    betat1(:) = betanornew(mayorage_elec+1,imun,1,ipauf,ieducf,elap,:)
                    betat2(:) = betanornew(mayorage_elec+1,imun,2,ipauf,ieducf,elap,:)
!!$                    write(*,*) 'beta1', betat1(:)
!!$                    write(*,*) 'beta2', betat2(:)                   
                 ELSE
                    betat1(:) = betanornew(mayorage_elec+1,imun,1,ipauf,ieducf,elap,:)
                 END IF
                 MMg = Nwealth
              END IF

              !We assign the probability of being audited, which depends on the audit regime ipau
              !It is used in mayorsolution0.f90 and mayorsolution.f90 to account for the probability of being audited the next term
              paudit_g = paudit(0,ipauf)
              IF (flag_audit_nterm == 1 .AND. nterms == Nnterms) paudit_g = paudit_nterm(0)

              CALL mayordecision(mayorage_elec,imun,ieducf,ability,ipauf,mwagef,is,msavingf,ifu,fundsf,ipr,zzprf,nterms,newpconsf,newqconsf,zzpuf,mnewsavingf,newstealf,utilf,vemaxpf)

!!$              if (myrank == 0) write(*,'(a,4f20.10)') 'frac_fine_g', frac_fine_c, newstealf, fine_fun(newstealf), paudit_g
              
              IF (mayorage_elec < Retage) THEN
                 DEALLOCATE(betat1,betat2)
              END IF
              IF (mayorage_elec < Retage-1 .AND. nterms < Nnterms .AND. (flag_no_run == 1 .OR. flag_norun_prob_policy == 1)) DEALLOCATE(betat3)

              !We save optimal stealing without measurement errors
              newstealf_nome = newstealf
              !We add a measurement error to stealing: A municipality can make a clerical/administrative mistake that will be detected by the auditors
              !We draw the measurement error from a normal distribution
!!$              meas_err_shock = (sigma_steal*rnormprob(i,j,tt)+mu_steal)*fundsf
!!$              IF (newstealf + meas_err_shock > 0d0) THEN
!!$                 newstealf = newstealf + meas_err_shock
!!$              ELSE
!!$                 newstealf = 0d0
!!$              END IF
              meas_err_shock = EXP(sigma_steal*rnormprob(i,j,tt)+mu_steal)
              newstealf = newstealf*meas_err_shock

!!$              if (myrank_t == 0 .AND. i <= 100 .AND. j <= 100) write(*,'(a,4i5,3f20.10)') 'meas errors', i, j, tt, nterms, newstealf_nome, newstealf, meas_err_shock
!!$              if (myrank_t == 0 .AND. i <= 100 .AND. j <= 100 .AND. nterms == 1) write(*,'(a,4i5,5f20.10)') 'ability, qcons', i, j, tt, nterms, mu_a, ability, newqconsf, zzpuf, newstealf
                               
              will_to_pay = 0d0
              ! We compute the value of the policy
              IF (flag_policy_evaluation == 1) THEN

                 ! We should compute the value function for the median resident only in the first period we consider
                 ! Since we compute the value function, it should include the effect of the policy for the rest of the median resident's life
                 ! THIS IS NOT TRUE BECAUSE IN tt == 1 ALL MAYORS ARE FIRST TERM AND THEY MAKE DIFFERNET DECISIONS
                 ! We use median age, education, and wealth conditional on the size of the municipality
                 med_edu_i = med_edu(imun)
                 med_ipast_i = med_ipast(imun)
                 med_age_g = med_age(imun)
                 med_wealth_g = med_wealth(imun)
                 pmwage_g = wage_vot(imun)
!!$                 pmwage_g = wage_pastmayor(med_age_g,imun,med_edu_i,med_ipast_i)
                 elapv = 1
                 IF (pappeal(1) < 0.5d0) elapv = 2

                 IF (ABS(delta-1d0) < small_no) THEN
                    ELNqcons = LOG(newqconsf/((popsizeg)**eta))
                 ELSE
                    ELNqcons = (newqconsf/((popsizeg)**eta))**(1d0-delta)/(1d0-delta)
                 END IF
                 ELNqcons_g =LOG(newqconsf/((popsizeg)**eta))

                 newsteal_g = newstealf
                 new_fr_stolen_g = newstealf/fundsf

                 !Find the point on the saving domain where we need to compute the approximated value function
                 DO isav = 1, Nwealth
                    IF (mnewsavingf <  wealth(isav)) THEN
                       ub_savi = isav
                       EXIT
                    END IF
                    ub_savi = Nwealth
                 END DO
                 
                 IF (nterms < Nnterms .AND. mayorage_elec < Retage-1) THEN
                    ALLOCATE(betat1(MMage_steal_q), betat2(MMage_steal_q))
                    MMg = MMage_steal_q
                    betat1(:) = betapmayor(med_age_g+1,imun,nterms+1,1,ieducf,med_ipast_i,ub_savi,ipauf,med_edu_i,elapv,:)
                    betat2(:) = betapmayor(med_age_g+1,imun,nterms+1,2,ieducf,med_ipast_i,ub_savi,ipauf,med_edu_i,elapv,:)
                    IF (flag_no_run == 1 .OR. flag_norun_prob_policy == 1) THEN
                       ALLOCATE(betat3(Nwealth))
                       betat3(:) = betapmayornor(med_age_g+1,imun,med_ipast_i,ipauf,med_edu_i,elapv,:)
                    END IF
                 ELSE IF (nterms >= Nnterms .OR. mayorage_elec >= Retage -1) THEN
                    ALLOCATE(betat1(Nwealth), betat2(Nwealth))
                    MMg = Nwealth
                    betat1(:) = betapmayornor(med_age_g+1,imun,med_ipast_i,ipauf,med_edu_i,elapv,:)
                    betat2(:) = betat1(:)
                 END IF

                 ! term of the current mayor; this is used in pastmayorsolution
                 itermi_g = nterms
                 agemayort = mayorage_elec

                 IF (flag_base_policy == 1) THEN

                    CALL notmayordecision(med_age_g,pmwage_g,med_wealth_g,ELNqcons)
                    ! We get the average value function for the municipality where the average is over simulations 
                    notmayor_value_t(i,j,tt) = vemaxoptpmg
                    ELNqcons_base_vec_t(i,j,tt) = ELNqcons

                    av_fr_stolen_t(i,tt) = av_fr_stolen_t(i,tt) + new_fr_stolen_g
                    terms_t(i,j,tt) = nterms
                    
                    IF (nterms == 1) THEN
                       fr_stolen_terms_t(1,i,j,tt) = new_fr_stolen_g
                    ELSE IF (nterms == 2) THEN
                       fr_stolen_terms_t(2,i,j,tt) = new_fr_stolen_g
                    ELSE IF (nterms == 3) THEN
                       fr_stolen_terms_t(3,i,j,tt) = new_fr_stolen_g
                    END IF

!!$                    IF (myrank_t >= 0) write(*,'(a,8i5,5f18.8)') 'notmayor_value_t a', i,j,tt,imun,nterms,med_age_g,med_edu_i,med_ipast_i,pmwage_g,med_wealth_g,ELNqcons_g,notmayor_value_t(i,j,tt),vemaxoptpmg
!!$                    IF (myrank_t == 0) write(*,*) 'betat', betat(:)

                 ELSE

                    IF (cost_policy == 1) THEN
                       IF (flag_prau_policy == 1 .OR. flag_optimal_prau_policy == 1 .OR. flag_term_optimal_prau_policy == 1 .OR. flag_norun_optimal_prau_policy == 1 .OR. flag_term_norun_optimal_prau_policy == 1) med_wealth_g = med_wealth_g - pc_cost_prau_t(imun)
                       IF (flag_wage_policy == 1) med_wealth_g = med_wealth_g - pc_cost_wage_t(imun)
                       IF (flag_optimal_wage_policy == 1 .OR. flag_wage_optimal_prau_policy == 1) THEN
                          med_wealth_g = med_wealth_g - pc_cost_prau_t(imun) - pc_cost_wage_t(imun)
                       END IF
                    END IF
                    
                    !We assign the citizens' value in the base case
                    notmayor_base_value = notmayor_value(i,j,tt)
                    ELNqcons_base = ELNqcons_base_vec(i,j,tt)

!!$596                 CALL notmayordecision(med_age_g,pmwage_g,med_wealth_g,ELNqcons_base)
!!$                    vemax_base = vemaxoptpmg
                    vemax_base = notmayor_base_value
                    
                    CALL notmayordecision(med_age_g,pmwage_g,med_wealth_g,ELNqcons)
                    vemax_policy = vemaxoptpmg

                    IF (vemax_base < vemax_policy) THEN
                       !Value function for the case with lower value, the value we want to achieve by lowering wealth
                       notmayor_low_value_g = vemax_base
                       !Value function for the other case
                       vemax_ini = vemax_policy
                       !Public consumption for the case with higher value, the public consumption we use when we compute the value function at lower values of wealth
                       ELNqcons_high_g = ELNqcons
                    ELSE
                       !Value function for the case with lower value, the value we want to achieve by lowering wealth
                       notmayor_low_value_g = vemax_policy
                       !Value function for the other case
                       vemax_ini = vemax_base
                       !Public consumption for the case with higher value, the public consumption we use when we compute the value function at lower values of wealth
                       ELNqcons_high_g = ELNqcons_base
                    END IF

!!$                       IF (flag_fafb == 1) write(*,'(A,2F20.10)') 'notmayor_base_value_g', notmayor_base_value_g, vemaxoptpmg
!!$                       IF (flag_fafb == 1) write(*,'(a,6i4,4f20.10)') 'notmayor_value_t a', i,j,med_age_g,imun,med_edu_i,med_ipast_i,pmwage_g,med_wealth_g,ELNqcons,ELNqcons_base

!!$                    write(*,'(A,2F20.10)') 'notmayor_base_value_g', notmayor_base_value_g, vemaxoptpmg
!!$                    IF (flag_fafb == 1) write(*,'(a,8i6)') 'notmayor_value_t dis', i,j,tt,imun,nterms,med_age_g,med_edu_i,med_ipast_i
!!$                    IF (flag_fafb == 1) write(*,'(a,7f30.20)') 'notmayor_value_t cont', pmwage_g,med_wealth_g,ELNqcons_base,ELNqcons,notmayor_base_value,vemax_base,vemax_ini
!!$                    IF (flag_fafb == 1) write(*,'(a,10f20.10)') 'betat', betat(:)
!!$                    IF (myrank_t >= 0) write(*,'(a,8i6)') 'notmayor_value_t dis', i,j,tt,imun,nterms,med_age_g,med_edu_i,med_ipast_i
!!$                    IF (myrank_t >= 0) write(*,'(a,7f30.20)') 'notmayor_value_t cont', pmwage_g,med_wealth_g,ELNqcons,ELNqcons_base,notmayor_base_low_g,vemax_base,vemax_policy
!!$                    IF (myrank_t >= 0) write(*,'(a,12i5,5f18.8)') 'mayor info',i,j,tt,mayorage_elec,imun,ieducf,ipauf,ipastf,is,ifu,ipr,nterms,ability,mwagef,msavingf,fundsf,zzprf

                    fr_wealth = 1d0
592                 a = med_wealth_g*fr_wealth
                    b = med_wealth_g
                    
                    fa = wealth_equivalent(a)
                    fb = wealth_equivalent(b)
                    
!!$                    IF (flag_fafb == 1) write(*,'(A,5F20.10)') 'fa, fb', a, fa, b, fb, fr_wealth
                    
                    IF(fa*fb <= 0.0d0) THEN
                       
                       !we look for the interior solution
                       wealth_star = zeroin(wealth_equivalent, a, b, 1d-10, 1d-10)
                       
                    ELSE
                       
                       IF (fb >= 0d0) THEN
                          fr_wealth = fr_wealth-0.2d0
                       ELSE
                          fr_wealth = fr_wealth+0.2d0
                       END IF
!!$                          flag_fafb = 1
!!$                          IF (myrank_t == 0) write(*,'(A,5F20.10)') 'WRONG STARTING POINTS', a, fa, b, fb, fr_wealth
!!$                          write(*,'(A,5F20.10)') 'WRONG STARTING POINTS', a, fa, b, fb, fr_wealth
                       GO TO 592
                       
                    END IF

                    IF (vemax_base < vemax_policy) THEN                    
                       will_to_pay = med_wealth_g - wealth_star
                    ELSE 
                       will_to_pay = wealth_star - med_wealth_g
                    END IF
                    av_will_to_pay_t(i,tt) = av_will_to_pay_t(i,tt) + will_to_pay
                    av_fr_stolen_t(i,tt) = av_fr_stolen_t(i,tt) + new_fr_stolen_g
                    terms_t(i,j,tt) = nterms
                    
                    IF (nterms == 1) THEN
                       fr_stolen_terms_t(1,i,j,tt) = new_fr_stolen_g
                    ELSE IF (nterms == 2) THEN
                       fr_stolen_terms_t(2,i,j,tt) = new_fr_stolen_g
                    ELSE IF (nterms == 3) THEN
                       fr_stolen_terms_t(3,i,j,tt) = new_fr_stolen_g
                    END IF
!!$                    if (myrank_t == 0) THEN
!!$                    IF (myrank_t ==0 .AND. i == 1 .AND. tt == 2) THEN
!!$                       write(*,'(a,3i4,7f14.8,6f18.9)') 'will_to_pay', tt, j, nterms, paudit_g, will_to_pay, med_wealth_g, wealth_star, vemax_base, vemax_policy, ELNqcons,newqconsf,zzpuf,mnewsavingf,newstealf,notmayor_base_value,ELNqcons_base
!!$                       write(*,'(a,15f12.7)') 'betat1', betat1
!!$                       write(*,'(a,15f12.7)') 'betat2', betat2
!!$                       !IF (nterms < Nnterms .AND. mayorage_elec < Retage-1) write(*,'(a,15f12.7)') 'betat3', betat3
!!$                       write(*,*) ''
!!$                    END if


!!$                    IF (flag_fafb == 1) write(*,'(a,3i5,5f20.10)') 'will_to_pay', i,j,tt,will_to_pay, med_wealth_g, wealth_star, notmayor_low_value_g, vemax_ini
!!$                    IF (myrank_t == 0 .AND. Nnterms == 1) write(*,'(a,3i5,5f20.10)') 'will_to_pay', i,j,tt, will_to_pay, new_fr_stolen_g, av_will_to_pay_t(i,tt), av_fr_stolen_t(i,tt), SUM(fr_stolen_terms_t(1,i,:,tt))
!!$                    IF (myrank_t == 0 .AND. Nnterms == 2) write(*,'(a,3i5,6f20.10)') 'will_to_pay', i,j,tt, will_to_pay, new_fr_stolen_g, av_will_to_pay_t(i,tt), av_fr_stolen_t(i,tt), SUM(fr_stolen_terms_t(1,i,:,tt)), SUM(fr_stolen_terms_t(2,i,:,tt)) 
!!$                    IF (myrank_t == 0 .AND. Nnterms == 3) write(*,'(a,3i5,7f20.10)') 'will_to_pay', i,j,tt, will_to_pay, new_fr_stolen_g, av_will_to_pay_t(i,tt), av_fr_stolen_t(i,tt), SUM(fr_stolen_terms_t(1,i,:,tt)), SUM(fr_stolen_terms_t(2,i,:,tt)), SUM(fr_stolen_terms_t(3,i,:,tt))
!!$                       IF (ABS(will_to_pay)>500d0) THEN
!!$                          flag_fafb = 1
!!$                          GO TO 596
!!$                       END IF
!!$
!!$                    flag_fafb = 0

                 END IF
                 
                 IF (mayorage_elec < Retage) THEN
                    DEALLOCATE(betat1,betat2)
                 END IF
                 IF (mayorage_elec < Retage-1 .AND. nterms < Nnterms .AND. (flag_no_run == 1 .OR. flag_norun_prob_policy == 1)) DEALLOCATE(betat3)
                 
              END IF
                               
!!$              IF (myrank_t == 0 .AND. newstealf < 0d0) write(*,'(a,6f16.10)') 'newsteal meas. errors after', newstealf, meas_err_shock, newstealf_nome, sigma_steal, rnormprob(i,j,tt), sigma_steal*rnormprob(i,j,tt)+mu_steal
!!$              IF (myrank_t == 0) write(*,'(a,6f16.10)') 'newsteal meas. errors after', newstealf, meas_err_shock, newstealf_nome, sigma_steal, rnormprob(i,j,tt), sigma_steal*rnormprob(i,j,tt)+mu_steal

              !If there is an infeasible solution we give the subsistence level to the mayor
              IF (utilf < check_bignumber .OR. ISNAN(mnewsavingf)) THEN 
                 infeasible = 1
                 WRITE(*,*) 'infeasible solution', newqconsf,mnewsavingf,newstealf,utilf,mayorage,imun,ieducf,ability,ist,auditf,mwagef,msavingf,fundsf,zzprf,yy,tt
                 newqconsf = subqcons
                 mnewsavingf = 0.0d0
                 newstealf = -999d0
                 utilf = subutil(INT(mayorage))
                 yy = -77777777777777d0
                 outelect = -777777
!!$                 flag_here = 890
!!$                 go to 890
              END IF

              IF(tt == 3) THEN
                 !audit
!!$                 uud(5) = udraw(i,j,1,5)
                 pauditt(1:2) = paudit(0:1,ipau)
                 IF (flag_audit_nterm == 1 .AND. nterms == Nnterms) pauditt(1:2) = paudit_nterm(0:1)
                 CALL rmultinomial(uud(5), 2, pauditt, aud3)
                 auditf3 = aud3-1
!!$                 IF(newstealf > l_steal_prob/norm) thief3 = 1 
                 IF(newstealf/fundsf > l_fr_stolen) thief3 = 1 
              END IF

              !We save the variables related to the new mayor
              simout(i,tt,12) = msavingf
              simout(i,tt,13) = DBLE(ieducf)
              simout(i,tt,14) = ability
              simout(i,tt,15) = DBLE(outelect)
              simout(i,tt,16) = mnewsavingf 
              simout(i,tt,17) = newstealf
              !newqconsf is in levels
              simout(i,tt,18) = newqconsf
              simout(i,tt,19) = yy
              simout(i,tt,20) = DBLE(nterms)
              simout(i,tt,21) = zzpuf
              simout(i,tt,22) = DBLE(ipastf)
              simout(i,tt,23) = DBLE(imun)
              simout(i,tt,24) = will_to_pay
              simout(i,tt,25) = newstealf_nome
              simout(i,tt,26) = DBLE(auditf3)
              simout(i,tt,27) = DBLE(mayorage_elec)
              simout(i,tt,28) = popsize(imun)
              simout(i,tt,29) = stealf
              IF (tt > 1) THEN
                 simout(i,tt,30) = rnorm(i,j,tt)
              ELSE
                 simout(i,tt,30) = 0d0
              END IF
              simout(i,tt,31) = newpconsf
              simout(i,tt,32) = incumb_age
              simout(i,tt,33) = utilf !frac_fine_c
              simout(i,tt,34) = vemaxpf !potential_wage
              simout(i,tt,35) = d_run
              simout(i,tt,36) = popsize(imun)**eta

              !This is only useful if we write out data
              IF (flag_simulate_data == 1 .OR. (flag_policy_evaluation == 1 .AND. nboot==1)) THEN
                 simoutt(i,j,tt,12) = msavingf
                 simoutt(i,j,tt,13) = DBLE(ieducf)
                 simoutt(i,j,tt,14) = ability
                 simoutt(i,j,tt,15) = DBLE(outelect)
                 simoutt(i,j,tt,16) = mnewsavingf 
                 simoutt(i,j,tt,17) = newstealf
                 !newqconsf is in levels
                 simoutt(i,j,tt,18) = newqconsf
                 simoutt(i,j,tt,19) = yy
                 simoutt(i,j,tt,20) = DBLE(nterms)
                 simoutt(i,j,tt,21) = zzpuf
                 simoutt(i,j,tt,22) = DBLE(ipastf)
                 simoutt(i,j,tt,23) = DBLE(imun)
                 ! Extra for new moment: dummy for steal >0 for discrete level of stealing
                 simoutt(i,j,tt,24) = will_to_pay
                 simoutt(i,j,tt,25) = newstealf_nome
                 simoutt(i,j,tt,26) = DBLE(auditf3)
                 simoutt(i,j,tt,27) = DBLE(mayorage_elec)
                 simoutt(i,j,tt,28) = popsize(imun)
                 simoutt(i,j,tt,29) = stealf
                 IF (tt > 1) THEN
                    simoutt(i,j,tt,30) = rnorm(i,j,tt)
                 ELSE
                    simoutt(i,j,tt,30) = 0d0
                 END IF
                 simoutt(i,j,tt,31) = newpconsf
                 simoutt(i,j,tt,32) = incumb_age
                 simoutt(i,j,tt,33) = utilf !frac_fine_c
                 simoutt(i,j,tt,34) = vemaxpf !potential_wage
                 simoutt(i,j,tt,35) = d_run
                 simoutt(i,j,tt,36) = popsize(imun)**eta
              END IF

           END DO
           
        END DO

        tt = 3        
        ! We compute the numerator and denominator for all the moments
        ii = 0
        
        ii = ii + 1
        ! Average fraction of funds stolen conditional on second term mayor
        IF (DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0)) > 0d0) THEN
           numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0))
           dent(ii) = dent(ii) + 1d0
        END IF
        
        ii = ii + 1
        ! Average fraction of funds stolen conditional on first term mayor
        IF (DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0)) > 0d0) THEN
           numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0))
           dent(ii) = dent(ii) + 1d0
        END IF
        
        !The probability of winning reelection conditional on running
        i_num = 0
        i_den = 0
        DO i = 1, Nactmun           
           IF (simout(i,tt,35) > 0.5d0) THEN
              i_den = i_den + 1
              IF (simout(i,tt,15) > 0.5d0) THEN
                 i_num = i_num + 1
              END IF
           END IF
        END DO
        ii = ii + 1
        IF (DBLE(i_den) > 0d0) THEN
           numt(ii) = numt(ii) + DBLE(i_num)/DBLE(i_den)
           dent(ii) = dent(ii) + 1d0
        END IF
        
        !We compute the fraction of incumbents that are high appeal at the time of the election (simout(:,2:3,8) > 1.5d0) in this simulation
        !We should use the decision for period 3, which corresponds to the decision to run at the end of the term 2000-2004 in the data
        i_num = 0
        i_den = 0
        DO i = 1, Nactmun           
           IF (simout(i,tt,35) > 0.5d0 .AND. simout(i,tt,7)> 0.5d0 .AND. simout(i,tt,4) > 0d0) THEN
              i_den = i_den + 1
              IF (simout(i,tt,15) > 0.5d0) THEN
                 i_num = i_num + 1
              END IF
           END IF
        END DO
        ii = ii + 1
        IF (DBLE(i_den) > 0d0) THEN
           numt(ii) = numt(ii) + DBLE(i_num)/DBLE(i_den)
           dent(ii) = dent(ii) + 1d0
        END IF
        
        ii = ii + 1
        ! Moment to identify the value of being in power:
        ! Fraction of mayors who did not steal and run for re-election; we should use the decision for period 3, which corresponds to the decision to run at the end of the term 2000-2004 in the data 
        thres_frac_stolen = l_fr_stolen
        ! We need to consider only mayors for which it was feasible to run (nterms < Nnterms, etc): simout(:,tt,35) > -0.5d0
        IF (DBLE(COUNT(simout(:,tt,35) > -0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,4) <= thres_frac_stolen .AND. simout(:,tt,19) > -99999999999999d0)) > 0d0) THEN
           numt(ii) = numt(ii) + DBLE(COUNT(simout(:,tt,35) > 0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,4) <= thres_frac_stolen .AND. simout(:,tt,19) > -99999999999999d0))/DBLE(COUNT(simout(:,tt,35) > -0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,4) <= thres_frac_stolen .AND. simout(:,tt,19) > -99999999999999d0))
           dent(ii) = dent(ii) + 1d0
        END IF
        
        ii = ii + 1
        !Probability of running for reelection conditional on stealing a positive amount of funds
        !We should use the decision for period 3, which corresponds to the decision to run at the end of the term 2000-2004 in the data
        !We need to consider only mayors for which it was feasible to run (nterms < Nnterms, etc): simout(:,tt,35) > -0.5d0
        IF (DBLE(COUNT(simout(:,tt,35) > -0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,4) > thres_frac_stolen .AND. simout(:,tt,19) > -99999999999999d0)) > 0d0) THEN
           numt(ii) = numt(ii) + DBLE(COUNT(simout(:,tt,35) > 0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,4) >  thres_frac_stolen .AND. simout(:,tt,19) > -99999999999999d0))/DBLE(COUNT(simout(:,tt,35) > -0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,4) > thres_frac_stolen .AND. simout(:,tt,19) > -99999999999999d0))
           dent(ii) = dent(ii) + 1d0
        END IF
        
        ! 75th perc. of qcons
        jj = 0
        DO i = 1, Nactmun
           ! Conditional on being feasible
           IF (simout(i,tt,19) > -99999999999999d0) THEN
              jj = jj + 1
              ! We record per-capita public consumpiton in the previous term
              qconst(jj) = simout(i,tt,5)/simout(i,tt,28)
           END IF
        END DO
        ptest = 75d0
        p75_qcons(1:1) = centile(qconst(1:jj), ptest)

	i_num = 0
        i_den = 0
        DO i = 1, Nactmun
!!$           IF (simout(i,tt,35) > 0.5d0 .AND. simout(i,tt,7)> 0.5d0 .AND. simout(i,tt,5) > p75_qcons(1)) THEN
           IF (simout(i,tt,35) > 0.5d0 .AND. simout(i,tt,5)/simout(i,tt,28) > p75_qcons(1)) THEN
              i_den = i_den + 1
              IF (simout(i,tt,15) > 0.5d0) THEN
                 i_num = i_num + 1
              END IF
           END IF
        END DO
        ii = ii + 1
        IF (DBLE(i_den) > 0d0) THEN
           numt(ii) = numt(ii) + DBLE(i_num)/DBLE(i_den)
           dent(ii) = dent(ii) + 1d0
        END IF
        
        ii = ii + 1
        ! Average fraction of funds stolen conditional on being in small municipality
        av_fr_stolen_small = SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,23))==1 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,23))==1 .AND. simout(:,2:2,19) > -99999999999999d0))

        ! Average fraction of funds stolen conditional on being in large municipality
        av_fr_stolen_large = SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,23))==3 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,23))==3 .AND. simout(:,2:2,19) > -99999999999999d0))

        numt(ii) = numt(ii) + av_fr_stolen_small - av_fr_stolen_large
        
        dent(ii) = dent(ii) + 1d0

!!$        write(*,'(a,2f20.10)') 'av_fr_stolen mun', av_fr_stolen_small, av_fr_stolen_large

        !The probability of winning reelection conditional on running, being audited and not stealing
        i_num = 0
        i_den = 0
        DO i = 1, Nactmun           
           IF (simout(i,tt,35) > 0.5d0 .AND. simout(i,tt,7)> 0.5d0 .AND. simout(i,tt,4) <= 0d0) THEN
              i_den = i_den + 1
              IF (simout(i,tt,15) > 0.5d0) THEN
                 i_num = i_num + 1
              END IF
           END IF
        END DO
        ii = ii + 1
        IF (DBLE(i_den) > 0d0) THEN
           numt(ii) = numt(ii) + DBLE(i_num)/DBLE(i_den)
           dent(ii) = dent(ii) + 1d0
        END IF

        ii = ii + 1
        ! Unconditional probablity of running (unconditional on being audited)
        IF (DBLE(COUNT(simout(:,tt,35) > -0.5d0 .AND. simout(:,tt,19) > -99999999999999d0)) > 0d0) THEN
           numt(ii) = numt(ii) + DBLE(COUNT(simout(:,tt,35) > 0.5d0 .AND. simout(:,tt,19) > -99999999999999d0))/DBLE(COUNT(simout(:,tt,35) > -0.5d0 .AND. simout(:,tt,19) > -99999999999999d0))
           dent(ii) = dent(ii) + 1d0
        END IF

        ii = ii + 1
        ! Moment to identify the parameter sigma_steal
        ! Std Dev of fraction stolen
        ! We use period 3 because we consider stealing and auditing for the previous term
        i_num = 0
        i_den = 0
        mean_lnstealing = 0d0
        mean2_lnstealing = 0d0
        DO i = 1, Nactmun           
           IF (simout(i,3,7) > 0.5d0) THEN
              i_den = i_den + 1
              mean_lnstealing = mean_lnstealing + simout(i,3,4)
              mean2_lnstealing = mean2_lnstealing + simout(i,3,4)**2d0
           END IF
!!$           IF (simout(i,3,8) > 0d0 .AND. simout(i,3,7) > 0.5d0) THEN
!!$              i_den = i_den + 1
!!$              mean_lnstealing = mean_lnstealing + LOG(simout(i,3,8))
!!$              mean2_lnstealing = mean2_lnstealing + LOG(simout(i,3,8))**2d0
!!$           END IF
!!$           IF (simout(i,3,7) > 0.5d0) THEN
!!$              i_num = i_num + 1
!!$              mean_stealing = mean_lnstealing + simout(i,3,8)
!!$              mean2_stealing = mean2_lnstealing + simout(i,3,8)**2d0
!!$           END IF
        END DO
        IF (i_den > 0) THEN
           mean_lnstealing = mean_lnstealing/DBLE(i_den)
           mean2_lnstealing = mean2_lnstealing/DBLE(i_den)           
           numt(ii) = numt(ii) + (mean2_lnstealing - mean_lnstealing**2d0)**0.5d0
           dent(ii) = dent(ii) + 1d0
!!$           if (myrank_t == 0) write(*,'(a,3f20.10)') 'mean and sd stealing', mean_lnstealing, mean2_lnstealing, (mean2_lnstealing - mean_lnstealing**2d0)**0.5d0
        END IF
!!$        IF (i_num > 0) THEN
!!$           mean_stealing = mean_stealing/DBLE(i_num)
!!$           mean2_stealing = mean2_stealing/DBLE(i_num)           
!!$           if (myrank_t == 0) write(*,'(a,3f20.10)') 'mean and sd stealing 223.1113', mean_stealing, mean2_stealing, (mean2_stealing - mean_stealing**2d0)**0.5d0
!!$        END IF
        
!!$        ii = ii + 1
!!$        ! Moment to identify the parameter sigma_steal
!!$        ! Prob(Steal<=0 & Auditf==1)/Prob(Audit==1) = fraction of mayor who stole
!!$        ! We use period 3 because we consider stealing and auditing for the previous term
!!$        IF (DBLE(COUNT(simout(:,3,4)>-1d0)) > 0d0) THEN
!!$           numt(ii) = numt(ii) + DBLE(COUNT(simout(:,3,4) <= 0d0 .AND. simout(:,3,7) > 0.5d0))/DBLE(COUNT(simout(:,3,4)>-1d0 .AND. simout(:,3,7) > 0.5d0 ))
!!$           dent(ii) = dent(ii) + 1d0
!!$        END IF
        
        ii = ii + 1
        ! Average per-capita public consumption        
        IF (DBLE(COUNT(simout(:,2:2,19) > -99999999999999d0)) > 0d0) THEN
           numt(ii) = numt(ii) + SUM(simout(:,2:2,18)/simout(:,2:2,28), simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(simout(:,2:2,19) > -99999999999999d0))
           dent(ii) = dent(ii) + 1d0
        END IF

        ii = ii + 1
        ! Moment to identify the fine parameter sigma_a (it is probably better to use the average of stealing squared).
        ! Prob(Steal>0 & Auditf==1)/Prob(Audit==1) = fraction of mayor who stole
        IF (DBLE(COUNT(simout(:,2:2,17)>-1d0 .AND. simout(:,2:2,19) > -99999999999999d0)) > 0d0) THEN
           numt(ii) = numt(ii) + DBLE(COUNT(simout(:,2:2,17)>thres_frac_stolen .AND. simout(:,2:2,19) > -99999999999999d0))/DBLE(COUNT(simout(:,2:2,17)>-1d0 .AND. simout(:,2:2,19) > -99999999999999d0))
           dent(ii) = dent(ii) + 1d0
        END IF
        
        ii = ii + 1
        ! Median stealing
        jj = 0
        DO i = 1, Nactmun
           ! Conditional on being feasible
           IF (simout(i,2,19) > -99999999999999d0) THEN
              jj = jj + 1
              ! We record stealing in the previous term
              fr_stolent(jj) = simout(i,2,17)/simout(i,2,9)
              !if (myrank_t == 0) write(*,*) 'new_steal', jj, newstealt(jj)
           END IF
        END DO
        ptest = 50d0
        p50_frstolen(1:1) = centile(fr_stolent(1:jj), ptest)
        numt(ii) = numt(ii) + p50_frstolen(1)
        dent(ii) = dent(ii) + 1d0
        
        ! Probablity of running conditional on being audited
        ii = ii + 1
        IF (DBLE(COUNT(simout(:,tt,35) > -0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,19) > -99999999999999d0)) > 0d0) THEN
           numt(ii) = numt(ii) + DBLE(COUNT(simout(:,tt,35) > 0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,19) > -99999999999999d0))/DBLE(COUNT(simout(:,tt,35) > -0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,19) > -99999999999999d0))
           dent(ii) = dent(ii) + 1d0
        END IF

        ! MOMENTS NOT USED IN ESTIMATION
        
!!$        ii = ii + 1
!!$        ! Average fraction of funds stolen conditional on first term mayor of high theta
!!$        numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) >= 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) >= 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of funds stolen conditional on first term mayor of high theta
!!$        numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) < 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) < 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of funds stolen conditional on first term mayor of high theta
!!$        numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,36) >= 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,36) >= 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of funds stolen conditional on first term mayor of high theta
!!$        numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,36) < 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,36) < 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of funds stolen conditional on first term mayor of high theta
!!$        numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) > 0d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) > 0d0 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of high type first term mayors
!!$        numt(ii) = numt(ii) + COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) >= 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of high type second term mayors
!!$        numt(ii) = numt(ii) + COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,36) >= 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ! Probability of winning reelection conditional on high type
!!$        i_num = 0
!!$        i_den = 0
!!$        DO i = 1, Nactmun           
!!$           IF (simout(i,tt,35) > 0.5d0 .AND. simout(i,tt,36) >= 1.5d0) THEN
!!$              i_den = i_den + 1
!!$              IF (simout(i,tt,15) > 0.5d0) THEN
!!$                 i_num = i_num + 1
!!$              END IF
!!$           END IF
!!$        END DO
!!$        ii = ii + 1
!!$        numt(ii) = numt(ii) + DBLE(i_num)/DBLE(i_den)
!!$        dent(ii) = dent(ii) + 1d0

        
!!$        tt = 3        
!!$        ! We compute the numerator and denominator for all the moments
!!$        ii = 0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of funds stolen conditional on second term mayor
!!$        numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of funds stolen conditional on first term mayor
!!$        numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        !This moment identifies the prob of a high type: The probability of winning reelection conditional on running
!!$        i_num = 0
!!$        i_den = 0
!!$        DO i = 1, Nactmun           
!!$           IF (simout(i,tt,35) > 0.5d0) THEN
!!$              i_den = i_den + 1
!!$              IF (simout(i,tt,15) > 0.5d0) THEN
!!$                 i_num = i_num + 1
!!$              END IF
!!$           END IF
!!$        END DO
!!$        ii = ii + 1
!!$        IF (i_den > 0) THEN
!!$           numt(ii) = numt(ii) + DBLE(i_num)/DBLE(i_den)
!!$           dent(ii) = dent(ii) + 1d0
!!$        END IF
!!$        
!!$        !This moment identifies the prob of a high type
!!$        !We compute the fraction of incumbents that are high appeal at the time of the election (simout(:,2:3,8) > 1.5d0) in this simulation
!!$        !We should use the decision for period 3, which corresponds to the decision to run at the end of the term 2000-2004 in the data
!!$        i_num = 0
!!$        i_den = 0
!!$        DO i = 1, Nactmun           
!!$           IF (simout(i,tt,35) > 0.5d0 .AND. simout(i,tt,7)> 0.5d0 .AND. simout(i,tt,4) > 0d0) THEN
!!$              i_den = i_den + 1
!!$              IF (simout(i,tt,15) > 0.5d0) THEN
!!$                 i_num = i_num + 1
!!$              END IF
!!$           END IF
!!$        END DO
!!$        ii = ii + 1
!!$        IF (i_den > 0) THEN
!!$           numt(ii) = numt(ii) + DBLE(i_num)/DBLE(i_den)
!!$           dent(ii) = dent(ii) + 1d0
!!$        END IF
!!$        
!!$        ii = ii + 1
!!$        ! Moment to identify the value of being in power:
!!$        ! Fraction of mayors who did not steal and run for re-election; we should use the decision for period 3, which corresponds to the decision to run at the end of the term 2000-2004 in the data 
!!$        thres_frac_stolen = l_fr_stolen
!!$        ! We need to consider only mayors for which it was feasible to run (nterms < Nnterms, etc): simout(:,tt,35) > -0.5d0
!!$        numt(ii) = numt(ii) + DBLE(COUNT(simout(:,tt,35) > 0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,4) <= thres_frac_stolen .AND. simout(:,tt,19) > -99999999999999d0))/DBLE(COUNT(simout(:,tt,35) > -0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,4) <= thres_frac_stolen .AND. simout(:,tt,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$           
!!$        ii = ii + 1
!!$        !Probability of running for reelection conditional on stealing a positive amount of funds
!!$        !We should use the decision for period 3, which corresponds to the decision to run at the end of the term 2000-2004 in the data
!!$        !We need to consider only mayors for which it was feasible to run (nterms < Nnterms, etc): simout(:,tt,35) > -0.5d0
!!$        numt(ii) = numt(ii) + DBLE(COUNT(simout(:,tt,35) > 0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,4) >  thres_frac_stolen .AND. simout(:,tt,19) > -99999999999999d0))/DBLE(COUNT(simout(:,tt,35) > -0.5d0 .AND. simout(:,tt,7) > 0.5d0 .AND. simout(:,tt,4) > thres_frac_stolen .AND. simout(:,tt,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Moment to identify the fine parameter sigma_a (it is probably better to use the average of stealing squared).
!!$        ! Prob(Steal>0 & Auditf==1)/Prob(Audit==1) = fraction of mayor who stole
!!$        numt(ii) = numt(ii) + DBLE(COUNT(simout(:,2:2,17)>l_fr_stolen .AND. INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0))/DBLE(COUNT(simout(:,2:2,17)>-1d0 .AND. INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0


        
!!$        ii = ii + 1
!!$        ! Average fraction of funds stolen conditional on first term mayor of high theta
!!$        numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) >= 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) >= 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of funds stolen conditional on first term mayor of high theta
!!$        numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) < 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) < 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of funds stolen conditional on first term mayor of high theta
!!$        numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,36) >= 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,36) >= 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of funds stolen conditional on first term mayor of high theta
!!$        numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,36) < 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,36) < 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of funds stolen conditional on first term mayor of high theta
!!$        numt(ii) = numt(ii) + SUM(simout(:,2:2,17)/simout(:,2:2,9), INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) > 0d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) > 0d0 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of high type first term mayors
!!$        numt(ii) = numt(ii) + COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,36) >= 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ii = ii + 1
!!$        ! Average fraction of high type second term mayors
!!$        numt(ii) = numt(ii) + COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,36) >= 1.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0
!!$
!!$        ! Probability of winning reelection conditional on high type
!!$        i_num = 0
!!$        i_den = 0
!!$        DO i = 1, Nactmun           
!!$           IF (simout(i,tt,35) > 0.5d0 .AND. simout(i,tt,36) >= 1.5d0) THEN
!!$              i_den = i_den + 1
!!$              IF (simout(i,tt,15) > 0.5d0) THEN
!!$                 i_num = i_num + 1
!!$              END IF
!!$           END IF
!!$        END DO
!!$        ii = ii + 1
!!$        IF (i_den > 0) THEN
!!$           numt(ii) = numt(ii) + DBLE(i_num)/DBLE(i_den)
!!$           dent(ii) = dent(ii) + 1d0
!!$        END IF

!!$        ! pth-percentile of stealing
!!$        jj = 0
!!$        DO i = 1, Nactmun
!!$           ! Conditional on running
!!$           IF (INT(simout(i,tt,35))==1) THEN
!!$              jj = jj + 1
!!$              ! We record stealing in the previous term
!!$              newstealt(jj) = simout(i,tt,4)
!!$              !if (myrank_t == 0) write(*,*) 'new_steal', jj, newstealt(jj)
!!$           END IF
!!$        END DO
!!$        ptest = 75d0
!!$        p75_newsteal(1:1) = centile(newstealt(1:jj), ptest)
!!$        !med_newsteal(1:1) = 0.1d0
!!$        p75_var = p75_newsteal(1)
!!$
!!$        ii = ii + 1
!!$        !Probability of running for reelection conditional on stealing more than the 75th percentile of funds
!!$        !We should use the decision for period 3, which corresponds to the decision to run at the end of the term 2000-2004 in the data
!!$        ! We need to consider only mayors for which it was feasible to run (nterms < Nnterms, etc): simout(:,tt,35) > -0.5d0
!!$        numt(ii) = numt(ii) + DBLE(COUNT(simout(:,tt,35) > 0.5d0 .AND. simout(:,tt,7) > 0.5 .AND. simout(:,tt,4) > p75_var .AND. simout(:,tt,19) > -99999999999999d0))/DBLE(COUNT(simout(:,tt,35) > -0.5d0 .AND. simout(:,tt,7) > 0.5 .AND. simout(:,tt,4) > p75_var .AND. simout(:,tt,19) > -99999999999999d0))
!!$        dent(ii) = dent(ii) + 1d0

        IF (DBLE(COUNT(INT(simout(:,2:2,20))==1)) > 0d0) THEN
           other_var_numt(1) = other_var_numt(1) + SUM(LOG(simout(:,2:2,14)), INT(simout(:,2:2,20))==1)/DBLE(COUNT(INT(simout(:,2:2,20))==1))
           other_var_dent(1) = other_var_dent(1) + 1d0
        END IF

        IF (DBLE(COUNT(INT(simout(:,2:2,20))==2)) > 0d0) THEN
           other_var_numt(2) = other_var_numt(2) + SUM(LOG(simout(:,2:2,14)), INT(simout(:,2:2,20))==2)/DBLE(COUNT(INT(simout(:,2:2,20))==2))
           other_var_dent(2) = other_var_dent(2) + 1d0
        END IF

        IF (DBLE(COUNT(INT(simout(:,2:2,20))==1)) > 0d0) THEN
           other_var_numt(3) = other_var_numt(3) + SUM(simout(:,2:2,14), INT(simout(:,2:2,20))==1)/DBLE(COUNT(INT(simout(:,2:2,20))==1))
           other_var_dent(3) = other_var_dent(3) + 1d0
        END IF

        IF (DBLE(COUNT(INT(simout(:,2:2,20))==2)) > 0d0) THEN
           other_var_numt(4) = other_var_numt(4) + SUM(simout(:,2:2,14), INT(simout(:,2:2,20))==2)/DBLE(COUNT(INT(simout(:,2:2,20))==2))
           other_var_dent(4) = other_var_dent(4) + 1d0
        END IF

        IF (DBLE(COUNT(simout(:,2:2,19) > -99999999999999d0)) > 0d0) THEN
           other_var_numt(5) = other_var_numt(5) + SUM(simout(:,2:2,18)/simout(:,2:2,36), simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(simout(:,2:2,19) > -99999999999999d0))
           other_var_dent(5) = other_var_dent(5) + 1d0
        END IF

        IF (DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0)) > 0d0) THEN
           other_var_numt(6) = other_var_numt(6) + SUM(simout(:,2:2,18)/simout(:,2:2,36), INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==1 .AND. simout(:,2:2,19) > -99999999999999d0))
           other_var_dent(6) = other_var_dent(6) + 1d0
        END IF
        
        IF (DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0)) > 0d0) THEN
           other_var_numt(7) = other_var_numt(7) + SUM(simout(:,2:2,18)/simout(:,2:2,36), INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(INT(simout(:,2:2,20))==2 .AND. simout(:,2:2,19) > -99999999999999d0))
           other_var_dent(7) = other_var_dent(7) + 1d0
        END IF
        
        IF (DBLE(COUNT(simout(:,2:2,35) > 0.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)) > 0d0) THEN
           other_var_numt(8) = other_var_numt(8) + SUM(simout(:,2:2,15), simout(:,2:2,35) > 0.5d0 .AND. simout(:,2:2,19) > -99999999999999d0)/DBLE(COUNT(simout(:,2:2,35) > 0.5d0 .AND. simout(:,2:2,19) > -99999999999999d0))
           other_var_dent(8) = other_var_dent(8) + 1d0
        END IF
        
     END IF

  END DO

  IF (flag_policy_evaluation == 1) THEN
          
     IF (flag_base_policy == 1) THEN

        IF (flag_optimal_auditprob == 1) THEN
           CALL MPI_ALLREDUCE(notmayor_value_t(1,1,1), notmayor_value(1,1,1), Nactmun*Nsim*Nsimterms, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)        
           CALL MPI_ALLREDUCE(ELNqcons_base_vec_t(1,1,1), ELNqcons_base_vec(1,1,1), Nactmun*Nsim*Nsimterms, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
        ELSE
           CALL MPI_ALLREDUCE(notmayor_value_t(1,1,1), notmayor_value(1,1,1), Nactmun*Nsim*Nsimterms, MPI_REAL8, MPI_SUM, NEW_COMM, ierr)        
           CALL MPI_ALLREDUCE(ELNqcons_base_vec_t(1,1,1), ELNqcons_base_vec(1,1,1), Nactmun*Nsim*Nsimterms, MPI_REAL8, MPI_SUM, NEW_COMM, ierr)
        END IF
        DEALLOCATE(notmayor_value_t)
        DEALLOCATE(ELNqcons_base_vec_t)
        
     END IF
     
     av_will_to_pay = 0d0
     av_fr_stolen = 0d0
     fr_stolen_terms = 0d0
     terms = 0
     IF (flag_optimal_auditprob == 1) THEN
        CALL MPI_ALLREDUCE(av_will_to_pay_t(1,1), av_will_to_pay(1,1), Nactmun*Nsimterms, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)     
        CALL MPI_ALLREDUCE(av_fr_stolen_t(1,1), av_fr_stolen(1,1), Nactmun*Nsimterms, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)     
        CALL MPI_ALLREDUCE(fr_stolen_terms_t(1,1,1,1), fr_stolen_terms(1,1,1,1), Nnterms*Nactmun*Nsim*Nsimterms, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
        CALL MPI_ALLREDUCE(terms_t(1,1,1), terms(1,1,1), Nactmun*Nsim*Nsimterms, MPI_INT, MPI_SUM, MPI_COMM_WORLD, ierr)
     ELSE        
        CALL MPI_ALLREDUCE(av_will_to_pay_t(1,1), av_will_to_pay(1,1), Nactmun*Nsimterms, MPI_REAL8, MPI_SUM, NEW_COMM, ierr)     
        CALL MPI_ALLREDUCE(av_fr_stolen_t(1,1), av_fr_stolen(1,1), Nactmun*Nsimterms, MPI_REAL8, MPI_SUM, NEW_COMM, ierr)     
        CALL MPI_ALLREDUCE(fr_stolen_terms_t(1,1,1,1), fr_stolen_terms(1,1,1,1), Nnterms*Nactmun*Nsim*Nsimterms, MPI_REAL8, MPI_SUM, NEW_COMM, ierr)
        CALL MPI_ALLREDUCE(terms_t(1,1,1), terms(1,1,1), Nactmun*Nsim*Nsimterms, MPI_INT, MPI_SUM, NEW_COMM, ierr)
     END IF     
     av_will_to_pay = av_will_to_pay/DBLE(Nsim)
     av_fr_stolen = av_fr_stolen/DBLE(Nsim)
     
!!$     DO i = 1, Nactmun
!!$        DO tt = 1, Nsimterms
!!$           IF (myrank_t == 0) write(*,'(a,2i6,2f20.10)') 'av_will_to_pay, av_fr_stolen', i, tt, av_will_to_pay(i,tt), av_fr_stolen(i,tt)
!!$           DO j = 1, Nnterms
!!$              IF (myrank_t == 0 .AND. COUNT(terms(i,:,tt) == j) > 0) write(*,'(a,2i6,1f20.10)') 'av_fr_stolen_terms', i, tt, SUM(fr_stolen_terms(j,i,:,tt), terms(i,:,tt) == j)/DBLE(COUNT(terms(i,:,tt) == j))
!!$           END DO
!!$        END DO
!!$     END DO

     IF (myrank_t == 0) write(*,'(a,1i4,2f20.10)') 'av_will_to_pay and av_fr_stolen over municipalities', nboot, SUM(av_will_to_pay(:,:))/DBLE(Nactmun*Nsimterms), SUM(av_fr_stolen(:,:))/DBLE(Nactmun*Nsimterms)     
     DO j = 1, Nnterms
        IF (myrank_t == 0 .AND. COUNT(terms(:,:,:) == j) > 0) write(*,'(a,2i4,1f20.10)') 'av_fr_stolen_terms over municipalities', nboot, j, SUM(fr_stolen_terms(j,:,:,:), terms(:,:,:) == j)/DBLE(COUNT(terms(:,:,:) == j))
     END DO
     
     av_wtp = SUM(av_will_to_pay(:,2))/DBLE(Nactmun)
     av_frst = SUM(av_fr_stolen(:,2))/DBLE(Nactmun)
     ! For terms 1 and 2 we consider period tt=2, for term 3 we have to use period tt=3 because in period 2 there cannot be a term 3
     DO j = 1, 2
        av_frst_terms(j) = SUM(fr_stolen_terms(j,:,:,2), terms(:,:,2) == j)/DBLE(COUNT(terms(:,:,2) == j))
     END DO
     IF (Nnterms == 3) av_frst_terms(j) = SUM(fr_stolen_terms(j,:,:,3), terms(:,:,3) == j)/DBLE(COUNT(terms(:,:,3) == j))
     
     IF (myrank_t == 0 .AND. Nnterms == 1) write(*,'(a,1i4,3f20.10)') 'av_will_to_pay and av_fr_stolen over municipalities for period 2', nboot, av_wtp, av_frst, av_frst_terms(1)
     IF (myrank_t == 0 .AND. Nnterms == 2) write(*,'(a,1i4,4f20.10)') 'av_will_to_pay and av_fr_stolen over municipalities for period 2', nboot, av_wtp, av_frst, av_frst_terms(1:2)
     IF (myrank_t == 0 .AND. Nnterms == 3) write(*,'(a,1i4,5f20.10)') 'av_will_to_pay and av_fr_stolen over municipalities for period 2', nboot, av_wtp, av_frst, av_frst_terms(1:3)
     
     av_wtp_resources = 100*av_wtp/95.770454d0
     years_to_death = 75d0-33.5d0
     av_wtp_income = 100d0*(av_wtp/years_to_death)/46.0146d0
     
     IF (myrank_t == 0) write(*,'(a,1i4,2f20.10)') 'av_will_to_pay reasources and income', nboot, av_wtp_resources, av_wtp_income
     
     IF (first_proc == 1) THEN
        iboot = nboot - Nboot_min + 1
        bootstrap_vect(policy_pos,iboot,1) = DBLE(nboot)
        bootstrap_vect(policy_pos,iboot,2) = av_wtp
        bootstrap_vect(policy_pos,iboot,3) = av_wtp_resources
        bootstrap_vect(policy_pos,iboot,4) = av_wtp_income
        bootstrap_vect(policy_pos,iboot,5) = av_frst
        bootstrap_vect(policy_pos,iboot,6) = av_frst_terms(1)       
        IF (Nnterms >= 2) bootstrap_vect(policy_pos,iboot,7) = av_frst_terms(2)
        IF (Nnterms == 3) bootstrap_vect(policy_pos,iboot,8) = av_frst_terms(3)     
        IF (myrank_t == 0) write(*,'(a,3i4,8f20.10)') 'bootstrap_vect', policy_pos, iboot, nboot, bootstrap_vect(policy_pos,iboot,:)
     END IF

!!$     IF (nboot==1 .AND. flag_base_policy == 0 .AND. flag_optimal_auditprob == 0) THEN
!!$
!!$        ALLOCATE(simouttt(Nactmun,Nsim,Nsimterms,Nvar))
!!$        simouttt = 0d0
!!$        
!!$        IF (flag_policy_evaluation == 0 .OR. (flag_policy_evaluation == 1 .AND. flag_optimal_auditprob == 1)) THEN
!!$           CALL MPI_ALLREDUCE(simoutt(1,1,1,1), simouttt(1,1,1,1), Nactmun*Nsim*Nsimterms*Nvar, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
!!$        ELSE
!!$           CALL MPI_ALLREDUCE(simoutt(1,1,1,1), simouttt(1,1,1,1), Nactmun*Nsim*Nsimterms*Nvar, MPI_REAL8, MPI_SUM, NEW_COMM, ierr)
!!$        END IF
!!$        
!!$        IF (myrank_t == 0) THEN
!!$           
!!$           ! POLICIES OUTPUT FILE
!!$           !INDIVIDUAL POLICIES
!!$           !Output file when we simulate the policy that increases the audit probability to 0.168 (Lula's policy)
!!$           IF (flag_prau_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_prau_Apr2024.dat", STATUS='REPLACE')
!!$           !Output file when we simulate the policy that sets the probability of releelection to zero if someone is caught stealing (CRA policy)
!!$           IF (flag_norun_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_norun_Apr2024.dat", STATUS='REPLACE')
!!$           !Output file when we simulate a policy in which mayors who stole in the past are audited with probability 1
!!$           !Output file when we simulate the policy that increases the term limit to 3
!!$           IF (flag_terms_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_terms_Apr2024.dat", STATUS='REPLACE')
!!$           !Output file when we simulate the policy that doubles wages
!!$           IF (flag_wage_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_wages_Apr2024.dat", STATUS='REPLACE')
!!$           
!!$           ! POLICY BUNDLES
!!$           !Output file when we simulate the policy that increases the term limit to 3 and sets the probability of releelection to zero if someone is caught stealing (3 TERMS PLUS CRA)
!!$           IF (flag_terms_norun_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_terms_norun_Apr2024.dat", STATUS='REPLACE')
!!$           !Output file when we simulate the three-term policy with higher probability of audit in the last term
!!$           IF (flag_terms_3th_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_terms_3termhigher_Apr2024.dat", STATUS='REPLACE')
!!$           !Output file when we simulate the three-term policy plus no run with higher probability of audit in the last term
!!$           IF (flag_terms_norun_3th_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_terms_norun_3termhigher_Apr2024.dat", STATUS='REPLACE')
!!$           !Output file when we simulate the policy that increases the audit probability to 0.168 (Lula's policy) jointly with the audited-if-caught policy 
!!$           !Output file when we simulate the no-run policy with higher probability of audit in the second term
!!$           IF (flag_norun_2th_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_norun_2termhigher_Apr2024.dat", STATUS='REPLACE')
!!$           !Output file when we simulate the no-run policy with higher probability of audit in the second term
!!$           IF (flag_terms_norun_wage_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_terms_norun_wages_Apr2024.dat", STATUS='REPLACE')
!!$           !Output file when we simulate the no-run policy with higher probability of audit in the second term
!!$           IF (flag_terms_norun_prauk_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_terms_norun_prauk_Apr2024.dat", STATUS='REPLACE')
!!$           
!!$           DO j = 1, Nsim
!!$              DO i = 1, Nactmun
!!$                 DO tt = 1, Nsimterms
!!$                    WRITE(28,*) simouttt(i,j,tt,:)
!!$                 END DO
!!$              END DO
!!$           ENDDO
!!$           CLOSE(28)
!!$
!!$           DEALLOCATE(simouttt)
!!$        
!!$        END IF
!!$     END IF

  END IF
  
  IF (flag_simulate_data == 1 .AND. flag_policy_evaluation == 0) THEN
     
     ALLOCATE(simouttt(Nactmun,Nsim,Nsimterms,Nvar))
     simouttt = 0d0
     
     IF (flag_policy_evaluation == 0 .OR. (flag_policy_evaluation == 1 .AND. flag_optimal_auditprob == 1)) THEN
        CALL MPI_ALLREDUCE(simoutt(1,1,1,1), simouttt(1,1,1,1), Nactmun*Nsim*Nsimterms*Nvar, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
     ELSE
        CALL MPI_ALLREDUCE(simoutt(1,1,1,1), simouttt(1,1,1,1), Nactmun*Nsim*Nsimterms*Nvar, MPI_REAL8, MPI_SUM, NEW_COMM, ierr)
     END IF

     IF (myrank_t == 0) THEN
     
        ! ESTIMATION OUTPUT FILE
        IF (flag_estimation == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_estimation_Apr2024.dat", STATUS='REPLACE')

        ! BASE CASE
        IF (flag_base_policy == 1 .AND. flag_no_selec_mun == 0 .AND. flag_no_selec_funds == 0 .AND. flag_no_selec_zzpr == 0 .AND. flag_no_selec_age == 0 .AND. flag_no_selec_educ == 0 .AND. flag_no_selec_ab == 0 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_base_Apr2024.dat", STATUS='REPLACE')

        ! POLICIES
        IF (flag_norun_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_norun_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate a policy in which mayors who stole in the past are audited with probability 1
        !Output file when we simulate the policy that increases the term limit to 3
        IF (flag_terms_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_terms_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate the policy that doubles wages
        IF (flag_wage_policy == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_wages_Apr2024.dat", STATUS='REPLACE')
        
        ! CONTROLLING FOR HETEROGENEITY
        !Output file when we simulate the policy that decreases the term limit to 1 
        IF (flag_1term_policy == 1 .AND. flag_no_selec_mun == 0 .AND. flag_no_selec_funds == 0 .AND. flag_no_selec_zzpr == 0 .AND. flag_no_selec_age == 0 .AND. flag_no_selec_educ == 0 .AND. flag_no_selec_ab == 0 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_1term_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate with no selection on municipality effects
        IF (flag_base_policy == 1 .AND. flag_no_selec_mun == 1 .AND. flag_no_selec_funds == 1 .AND. flag_no_selec_zzpr == 1 .AND. flag_no_selec_age == 0 .AND. flag_no_selec_educ == 0 .AND. flag_no_selec_ab == 0 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_no_sel_mun_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate with no selection on ability
        IF (flag_base_policy == 1 .AND. flag_no_selec_mun == 0 .AND. flag_no_selec_funds == 0 .AND. flag_no_selec_zzpr == 0 .AND. flag_no_selec_age == 0 .AND. flag_no_selec_educ == 0 .AND. flag_no_selec_ab == 1 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_no_sel_ab_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate with no selection on educ
        IF (flag_base_policy == 1 .AND. flag_no_selec_mun == 0 .AND. flag_no_selec_funds == 0 .AND. flag_no_selec_zzpr == 0 .AND. flag_no_selec_age == 0 .AND. flag_no_selec_educ == 1 .AND. flag_no_selec_ab == 0 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_no_sel_edu_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate with no selection on educ
        IF (flag_base_policy == 1 .AND. flag_no_selec_mun == 0 .AND. flag_no_selec_funds == 0 .AND. flag_no_selec_zzpr == 0 .AND. flag_no_selec_age == 1 .AND. flag_no_selec_educ == 0 .AND. flag_no_selec_ab == 0 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_no_sel_age_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate with no selection on human capital
        IF (flag_base_policy == 1 .AND. flag_no_selec_mun == 0 .AND. flag_no_selec_funds == 0 .AND. flag_no_selec_zzpr == 0 .AND. flag_no_selec_age == 0 .AND. flag_no_selec_educ == 1 .AND. flag_no_selec_ab == 1 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_no_sel_hc_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate with no selection on human capital and age
        IF (flag_base_policy == 1 .AND. flag_no_selec_mun == 0 .AND. flag_no_selec_funds == 0 .AND. flag_no_selec_zzpr == 0 .AND. flag_no_selec_age == 1 .AND. flag_no_selec_educ == 1 .AND. flag_no_selec_ab == 1 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_no_sel_hcage_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate with no selection on savings
        IF (flag_base_policy == 1 .AND. flag_no_selec_mun == 0 .AND. flag_no_selec_funds == 0 .AND. flag_no_selec_zzpr == 0 .AND. flag_no_selec_age == 0 .AND. flag_no_selec_educ == 0 .AND. flag_no_selec_ab == 0 .AND. flag_no_selec_sav == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_no_sel_sav_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate the policy that imposes a 1 term limit and we only remove selection on municipality effects
        IF (flag_1term_policy == 1 .AND. flag_no_selec_mun == 1 .AND. flag_no_selec_funds == 1 .AND. flag_no_selec_zzpr == 1 .AND. flag_no_selec_age == 0 .AND. flag_no_selec_educ == 0 .AND. flag_no_selec_ab == 0 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_1term_no_sel_mun_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate the policy that imposes a 1 term limit and we only remove selection on human capital
        IF (flag_1term_policy == 1 .AND. flag_no_selec_mun == 0 .AND. flag_no_selec_funds == 0 .AND. flag_no_selec_zzpr == 0 .AND. flag_no_selec_age == 0 .AND. flag_no_selec_educ == 1 .AND. flag_no_selec_ab == 1 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_1term_no_sel_hc_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate the policy that imposes a 1 term limit and we only remove selection on human capital and age
        IF (flag_1term_policy == 1 .AND. flag_no_selec_mun == 0 .AND. flag_no_selec_funds == 0 .AND. flag_no_selec_zzpr == 0 .AND. flag_no_selec_age == 1 .AND. flag_no_selec_educ == 1 .AND. flag_no_selec_ab == 1 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_1term_no_sel_hcage_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate the policy that imposes a 1 term limit and we only remove selection on savings
        IF (flag_1term_policy == 1 .AND. flag_no_selec_mun == 0 .AND. flag_no_selec_funds == 0 .AND. flag_no_selec_zzpr == 0 .AND. flag_no_selec_age == 0 .AND. flag_no_selec_educ == 0 .AND. flag_no_selec_ab == 0 .AND. flag_no_selec_sav == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_1term_no_sel_sav_Apr2024.dat", STATUS='REPLACE')
        !Output file when we simulate the policy that imposes a 1 term limit and we only remove selection on municipality effects, human capital, and age
        IF (flag_1term_policy == 1 .AND. flag_no_selec_mun == 1 .AND. flag_no_selec_funds == 1 .AND. flag_no_selec_zzpr == 1 .AND. flag_no_selec_age == 1 .AND. flag_no_selec_educ == 1 .AND. flag_no_selec_ab == 1 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_1term_no_sel_obs_unobs_Apr2024.dat", STATUS='REPLACE')
        IF (flag_base_policy == 1 .AND. flag_no_selec_mun == 1 .AND. flag_no_selec_funds == 1 .AND. flag_no_selec_zzpr == 1 .AND. flag_no_selec_age == 1 .AND. flag_no_selec_educ == 1 .AND. flag_no_selec_ab == 1 .AND. flag_no_selec_sav == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_no_sel_obs_unobs_sav_Apr2024.dat", STATUS='REPLACE')
        IF (flag_base_policy == 1 .AND. flag_no_selec_mun == 1 .AND. flag_no_selec_funds == 1 .AND. flag_no_selec_zzpr == 1 .AND. flag_no_selec_age == 1 .AND. flag_no_selec_educ == 1 .AND. flag_no_selec_ab == 1 .AND. flag_no_selec_sav == 0) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_no_sel_obs_unobs_Apr2024.dat", STATUS='REPLACE')
        IF (flag_1term_policy == 1 .AND. flag_no_selec_age == 1 .AND. flag_no_selec_educ == 1 .AND. flag_no_selec_ab == 1 .AND. flag_no_selec_sav == 1 .AND. flag_no_selec_funds == 1 .AND. flag_no_selec_zzpr == 1 .AND. flag_no_selec_mun == 1) OPEN(unit=28, FILE="~/MyPapers/Mayors/Estimation/Stata/simdata_1term_no_sel_Apr2024.dat", STATUS='REPLACE')
        
        DO j = 1, Nsim
           DO i = 1, Nactmun
              DO tt = 1, Nsimterms
                 WRITE(28,*) simouttt(i,j,tt,:)
              END DO
           END DO
        ENDDO
        CLOSE(28)
        
        DEALLOCATE(simouttt)
        
     END IF
  END IF
     
  IF (flag_policy_evaluation == 0 .OR. (flag_policy_evaluation == 1 .AND. flag_optimal_auditprob == 1)) THEN
     CALL MPI_BARRIER(MPI_COMM_WORLD, ierr)
  ELSE
     CALL MPI_BARRIER(NEW_COMM, ierr)
  END IF
  
!!$     DEALLOCATE(notmayor_value, ELNqcons_base_vec)     
     
  IF(myrank_t == 0) write(*,*) ' ------------------------------------------------------------------'
  IF(myrank_t == 0) write(*,*) '  Compute the simulated moments for the estimation                 '
  IF(myrank_t == 0) write(*,*) ' ------------------------------------------------------------------'

  IF (flag_policy_evaluation == 0 .OR. (flag_policy_evaluation == 1 .AND. flag_optimal_auditprob == 1)) THEN

     CALL MPI_BARRIER(MPI_COMM_WORLD, ierr)
     CALL MPI_ALLREDUCE(numt(1), num(1), Nmomtemp, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
     CALL MPI_ALLREDUCE(dent(1), den(1), Nmomtemp, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
     CALL MPI_ALLREDUCE(other_var_numt(1), other_var_num(1), N_other_var, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)
     CALL MPI_ALLREDUCE(other_var_dent(1), other_var_den(1), N_other_var, MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, ierr)

  ELSE

     CALL MPI_BARRIER(NEW_COMM, ierr)
     CALL MPI_ALLREDUCE(numt(1), num(1), Nmomtemp, MPI_REAL8, MPI_SUM, NEW_COMM, ierr)
     CALL MPI_ALLREDUCE(dent(1), den(1), Nmomtemp, MPI_REAL8, MPI_SUM, NEW_COMM, ierr)
     CALL MPI_ALLREDUCE(other_var_numt(1), other_var_num(1), N_other_var, MPI_REAL8, MPI_SUM, NEW_COMM, ierr)
     CALL MPI_ALLREDUCE(other_var_dent(1), other_var_den(1), N_other_var, MPI_REAL8, MPI_SUM, NEW_COMM, ierr)

  END IF

     
  ! Print moments used in the estimation
  DO i = 1, nmom
     
     IF (ABS(den(i)) > 0d0) THEN
        simmom(i,1) = num(i) / den(i)
     ELSE
        simmom(i,1) = 0d0
        IF (myrank_t == 0) write(*,'(a,1i4,3f20.10)') 'den = 0d0', i, simmom(i,1), num(i), den(i)
     END IF

     IF (ISNAN(simmom(i,1))) THEN
        simmom(i,1) = 100d0
        IF (myrank_t == 0) write(*,'(a,1i4,3f20.10)') 'simmom NaN', i, simmom(i,1), num(i), den(i)
     END IF
     
  END DO

  IF (nmomtemp > nmom) THEN
     ! Print moments not used in the estimation
     DO i = nmom+1, nmomtemp
        
        IF (ABS(den(i)) > 0d0) THEN
           simmom_p(i,1) = num(i) / den(i)
        ELSE
           simmom_p(i,1) = 0d0
           IF (myrank_t == 0) write(*,'(a,1i4,3f20.10)') 'den_p = 0d0', i, simmom_p(i,1), num(i), den(i)
        END IF
        
        IF (ISNAN(simmom_p(i,1))) THEN
           simmom_p(i,1) = 100d0
           IF (myrank_t == 0) write(*,'(a,1i4,3f20.10)') 'simmom_p NaN', i, simmom_p(i,1), num(i), den(i)
        END IF

     END DO
  END IF
  
  DO i = 1, N_other_var
     IF (ABS(other_var_den(i)) > 0d0) THEN
        other_var(i) = other_var_num(i) / other_var_den(i)
     ELSE
        other_var(i) = 0d0
     END IF
  END DO
  
  transform = 0d0  
  DO i = 1, nmom
     IF (flag_sim == 0) THEN
        datamom(i,1) = data(i,1)*norm_factor(i)
        simmom(i,1) = simmom(i,1)*norm_factor(i)
     ELSE
        !We use the moments from the simulated data
        datamom(i,1) = simdatamom(i,1)*norm_factor(i)
     END IF
     transform(i,i)= norm_factor(i)
  END DO
  IF (nmomtemp > nmom) THEN
     ! Print moments not used in the estimation
     DO i = nmom+1, nmomdata
        IF (flag_sim == 0) THEN
           datamom(i,1) = data(i,1)*norm_factor(i)
        ELSE
           !We use the moments from the simulated data
           datamom(i,1) = simdatamom(i,1)*norm_factor(i)
        END IF
     END DO
  END IF
  
  sim_var_names = [character(len = 100) :: 'Average Stealing in Second Term, SIMULATION', &
       'Average Stealing in First Term, SIMULATION', &
       'Probability Winning Reelection, SIMULATION', &
       'Probability Winning Reelection Conditional on Stealing, SIMULATION', &
       'Probability of Running Conditional on Not Stealing, SIMULATION', &
       'Probability of Running Conditional on Stealing, SIMULATION', &
       'Probability Winning Reelection Conditional on Qcons > 75th perc., SIMULATION', &
       'Difference in Average Stealing, Small Vs Large Mun, SIMULATION', &
       'Probability Winning Reelection Conditional on not Stealing, SIMULATION', &
       'Probability of Running, SIMULATION', &
       'Fraction of Mayor Not Stealing, SIMULATION', &            
       'Per-capita public consumption, SIMULATION', &
       'Median Stealing , SIMULATION', &
       'Probability of Running Conditional on Being Audited, SIMULATION', &
       'Average Stealing in First Term by High Type Mayors, SIMULATION', &
       'Average Stealing in First Term by Low Type Mayors, SIMULATION', &
       'Average Stealing in Second Term by High Type Mayors, SIMULATION']

  data_var_names = [character(len = 100) :: 'Average Stealing in Second Term, DATA', &
       'Average Stealing in First Term, DATA', &
       'Probability Winning Reelection, DATA', &
       'Probability Winning Reelection Conditional on Stealing More Than the Median, DATA', &
       'Probability of Running Conditional on Not Stealing, DATA', &
       'Probability of Running Conditional on Stealing, DATA', &
       'Probability Winning Reelection Conditional on Qcons > 75th perc., DATA', &
       'Difference in Average Stealing, Small Vs Large Mun, DATA', &
       'Probability Winning Reelection Conditional on not Stealing, DATA', &
       'Probability of Running, DATA', &
       'Fraction of Mayor Not Stealing, DATA', &     
       'Per-capita public consumption, DATA', &
       'Median Stealing , DATA', &
       'Probability of Running Conditional on Being Audited, DATA', &
       'Average Stealing in First Term by High Appeal Mayors, DATA', &
       'Average Stealing in First Term by Low Appeal Mayors, DATA', &
       'Average Stealing in Second Term by High Appeal Mayors, DATA']

  IF (myrank_t == 0) THEN
     i = 0
     DO ii = 1, nmomtemp
        i = i + 1
        write(*,*) ''
        IF (i <= nmom) THEN
           write(*,*) trim(sim_var_names(i)), simmom(i,1)
        ELSE
           write(*,*) trim(sim_var_names(i)), simmom_p(i,1)
        END IF
        IF (i <= nmomdata) write(*,*) trim(data_var_names(i)), datamom(i,1)
     END DO
     write(*,*) ''
     write(*,*) 'Average Log Ability 1st Term, SIMULATION', other_var(1)
     write(*,*) ''
     write(*,*) 'Average Log Ability 2nd Term, SIMULATION', other_var(2)
     write(*,*) ''
     write(*,*) 'Difference in Log Average Ability between 2nd and 1st Term, SIMULATION', other_var(2)-other_var(1)
     write(*,*) ''
     write(*,*) 'Average Ability 1st Term, SIMULATION', other_var(3)
     write(*,*) ''
     write(*,*) 'Average Ability 2nd Term, SIMULATION', other_var(4)
     write(*,*) ''
     write(*,*) 'Difference in Average Ability between 2nd and 1st Term, SIMULATION', other_var(4)-other_var(3)
     write(*,*) ''
     write(*,*) 'Per-capita public consumption, SIMULATION', other_var(5)
     write(*,*) ''
     write(*,*) 'Per-capita public consumption 1st term, SIMULATION', other_var(6)     
     write(*,*) ''
     write(*,*) 'Per-capita public consumption 2nd term, SIMULATION', other_var(7)
     write(*,*) ''
     write(*,*) 'Fraction Elected, SIMULATION', other_var(8)
     write(*,*) ''
  END IF
  
  DEALLOCATE(simout)
  flag_simulation = 0
  
  CALL cpu_time(ts2)
  IF (myrank_t == 0) write(*,*) 'time simulation', ts2-ts1

END SUBROUTINE simulatemun
   
FUNCTION wealth_equivalent(wealthin)
  USE commonvars
  IMPLICIT NONE
  REAL(8), INTENT(IN) :: wealthin
  REAL(8)             :: wealth_equivalent

  EXTERNAL            :: notmayordecision

  CALL notmayordecision(med_age_g,pmwage_g,wealthin,ELNqcons_high_g)

  wealth_equivalent = vemaxoptpmg - notmayor_low_value_g

!!$  IF (myrank_t == 0) write(*,'(a,4f20.10)') 'wealth_equivalent', wealth_equivalent, wealthin, vemaxoptpmg, notmayor_low_value_g
!!$  IF (flag_fafb == 1) write(*,'(a,4f20.10)') 'wealth_equivalent', wealth_equivalent, wealthin, vemaxoptpmg, notmayor_low_value_g

END FUNCTION wealth_equivalent
